Controle Automático de Aeronaves

Mecânica do Vôo e Controle de Aeronaves

•
UNESP

Alexandre Masson Vicente
25/03/2022
E aí, curtiu este material?
Ajude a incentivar outros estudantes a melhorar o conteúdo
Gostou desse material? Compartilhe! 🧡
Mecânica do Vôo e Controle de Aeronaves

11 Materiais compartilhados
Baixe o app para aproveitar ainda mais
Leia os materiais offline, sem usar a internet. Além de vários outros recursos!
Prévia do material em texto
MASSON
CONTROLE AUTOMÁTICO
DE
AERONAVES
Uberlândia
Alexandre Masson Vicente
2017
UNIVERSIDADE FEDERAL DE UBERLÂNDIA
APOSTILA
MECÂNICA DO VÔO
Controle Automático de Aeronaves
Autor:
Alexandre MASSON
Supervisor:
Dr. Leonardo SANCHES
2017
Resumo
Este documento contém estudos de caso em controle automático de aero-
naves. Ele cobre, inicialmente, o problema de busca de regiões de operação
em cima das quais o sistema de controle de uma aeronave rı́gida de asa fixa
deverá operar e, posteriormente, a linearização de suas equações de movi-
mento em torno de tais pontos. Serão extraı́dos deste processo duas classes
de modelos: funções de transferência e variáveis de estado. Em seguida,
será trabalhado o projeto de compensadores baseados na Teoria de Controle
Clássico para o veı́culo modelado dentro de suas dinâmicas latero-direcional
e longitudinal. Na sequência, para um quadricóptero, são desenvolvidos
controladores baseados na Teoria de Controle Moderno tanto para o caso
de regulação quanto para o rastreamento, explorando conceitos básicos do
assunto ao longo da resolução. E por fim, são pontuados materiais com-
plementares a esta apostila que podem ser combinados com a mesma para
potencializar o trabalho de profissionais do ensino de controle.
Sumário
1 Introdução 2
2 Modelagem Matemática da Aeronave 4
2.1 Sistemas de Referência . . . . . . . . . . . . . . . . . . . . . . . 4
2.2 Dinâmica do Corpo Rı́gido . . . . . . . . . . . . . . . . . . . . . 5
2.3 Dinâmica de Voo . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.4 Trimagem e Linearização do Modelo . . . . . . . . . . . . . . . . 11
2.5 Funções de Transferência de Regimes de Voo: O Caso da Curva
Coordenada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3 Projeto de Compensadores de Voo 19
3.1 Controle do voo latero-direcional . . . . . . . . . . . . . . . . . . 19
3.1.1 Rolagem . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.1.2 Ângulo de Derrapagem . . . . . . . . . . . . . . . . . . . 32
3.2 Controle do voo longitudinal . . . . . . . . . . . . . . . . . . . . 44
3.2.1 Rudimentos de Resposta em Frequência . . . . . . . . . . 44
3.2.2 Redes proporcional, em avanço de fase e em atraso de fase 61
3.2.2.1 Compensador proporcional . . . . . . . . . . . 62
3.2.2.2 Rede em avanço . . . . . . . . . . . . . . . . . 62
1
3.2.2.3 Rede em atraso . . . . . . . . . . . . . . . . . 64
3.2.3 Altitude e velocidade . . . . . . . . . . . . . . . . . . . . 72
3.2.3.1 Altitude . . . . . . . . . . . . . . . . . . . . . 74
3.2.3.2 Velocidade . . . . . . . . . . . . . . . . . . . . 82
4 Prelúdio ao Controle Moderno: o problema do quadricóptero 87
4.1 Modelo em espaço de estados do quadricóptero . . . . . . . . . . 88
4.2 A regulação de estados . . . . . . . . . . . . . . . . . . . . . . . 92
4.3 O rastreamento de estados . . . . . . . . . . . . . . . . . . . . . 111
5 Materiais complementares 130
5.1 Recursos didáticos auxiliares ao conteúdo da apostila . . . . . . . 130
5.2 Sugestões de literatura voltada ao ensino de controle automático . 131
5.3 Recomendações de leitura para aspirantes à área de Controle de
Sistemas Dinâmicos . . . . . . . . . . . . . . . . . . . . . . . . . 131
6 Conclusão 132
Bibliografia 132
1 Introdução
A navegação e guiagem de aeronaves é um problema de Engenharia Aeronáutica
encontrado em aplicações conhecidas nas aviações militar, comercial e para fins
de pesquisa: agricultura, mapeamento de florestas, reconhecimento e delimitação
de fronteiras territoriais etc. E sua execução é possı́vel através da integração de
vários sub-sistemas que caracterizam a Aviônica de uma aeronave, tais como:
• drivers de amplificação de sinais, sejam eles hidráulicos, pneumáticos ou de
eletrônica de potência);
• sensores de navegação inercial e filtros atenuadores de ruı́dos;
• superfı́cies provedoras de forças e torques de atuação;
• eletrônica embarcada de controle capaz de processar a saı́da dos sensores e
os comandos do piloto em sinais que entrarão nos atuadores.
2
Figura 1: Exemplo de um sistema de controle de vôo.
Este trabalho tem por objetivo explorar as leis de controle que são implemen-
tadas no último sub-sistema dos tópicos listados acima. A construção delas se dará
pela aplicação das teorias de Controle Clássico e Controle Moderno a problemas
de estabilização e regulação de aeronaves. Nestes, serão explorados conceitos fun-
damentais de Análise de Sistemas Dinâmicos como estabilidade, resposta tran-
sitória e resposta em regime permanente.
Além disto, será ilustrada a conexão entre o modelo matemático da aeronave e
os aspectos de análise com rudimentos da Teoria da Perturbação. E, por fim, al-
guns pontos práticos do projeto de técnicas controle serão elucidados nos estudos
de caso, por exemplo rejeição de distúrbios e/ou ruı́dos, saturação dos atuadores
etc.
Para cada resolução de um caso, serão obtidos resultados oriundos de simulações
numéricas, em ambiente MATLAB, com o intuito de verificar as predições da teo-
ria. Contudo, é pertinente salientar que o projeto de Controle Automático de
Aeronaves envolve todas as áreas do conhecimento mencionadas aqui em todas
3
as fases de seu desenvolvimento, o que ficará mais evidente no decorrer do tra-
balho.
2 Modelagem Matemática da Aeronave
Inicialmente, o foco da modelagem se dará em torno de uma aeronave de asa
fixa como a da Figura 4.
Visto que a massa de uma aeronave desse tipo é pequena, se comparada à
massa de objetos astronômicos como planetas, por exemplo, e que sua velocidade
é baixa, quando sua magnitude é confrontada com o tamanho da velocidade da luz,
efeitos relativı́sticos no seu movimento podem ser desprezados. Além disso, dadas
as dimensões da aeronave, os valores de energia e módulo do momento linear a
ela associados durante o movimento serão sempre altos o bastante para ignorar
efeitos de propagação ondulatória da matéria que a compõe durante sua trajetória.
Portanto, é pertinente obter o seu modelo matemático através da Mecânica Newto-
niana. Ou seja, o avião é um emaranhado contı́nuo de corpúsculos e as descrições
do movimento advindas da análise de observadores localizados em referenciais
galileanos serão equivalentes; ademais, a gravidade poderá ser tratada através de
um modelo de força externa.
Dentro da validade da Mecânica Clássica para a abordagem do problema,
serão assumidas, ainda, duas hipóteses: a aeronave, para todos os efeitos, é um
corpo rı́gido e a Primeira Lei de Newton vale para observações feitas da Terra,
dois pontos assumidos que permitirão a utilização dos Princı́pios de Newton-
Euler, da Dinâmica do Corpo Rı́gido, para a obtenção das equações diferenciais
do veı́culo. Hipóteses adicionais, juntamente com suas consequências lógicas,
serão exploradas no decorrer das sessões que seguem.
2.1 Sistemas de Referência
Para avaliar de forma satisfatória o movimento do veı́culo, devemos escolher
um local no qual seja possı́vel fazer medições de deslocamento e o orientação
espaciais do veı́culo. E a superfı́cie da Terra é um lugar razoável para se efetuar
essas observações.
Negligenciando efeitos de aceleração do planeta e de rotação em torno do
seu centro, pode-se assumir essa superfı́cie como um Referencial Inercial. Mas
como alguns sensores se encontram embarcados na aeronave, inevitavelmente, no
decorrer da formulação, será necessário representá-la, também, como uma trı́ade
4
de eixos centrada no seu baricentro, que será aqui denotada por Referencial do
Corpo.
Figura 2: Aeronave e solo.
Os sistemas inercial e do corpo podem ser denotados por suas respectivas
bases de vetores unitários, F e F ′ , definidas a seguir:
F = {ê1, ê2, ê3} (1)
F ′ = {ê′1, ê′2, ê′3} (2)
2.2 Dinâmica do Corpo Rı́gido
Uma vez definidos os referenciais em relação aos quais as quantidadescinemáticas
serão medidas e decompostas (em 3 direções), para avaliar o movimento da aero-
nave rı́gida, serão enunciados os Princı́pios de Newton-Euler, válidos no Referen-
cial Inercial:
F = Ṗ (3)
onde P é o momento linear do veı́culo e F é a resultante das forças externas
nele aplicadas, ambos medidos no Referencial Inercial.
N = L̇ (4)
onde L é o momento angular do veı́culo eN é a resultante dos torques exter-
nos nele aplicados, ambos medidos no Referencial Inercial.
5
Expressando a dinâmica de translação, representada por 3, em função da veloci-
dade de translação do veı́culo, tem-se que:
F = mv̇ = m(v̇v̂ + ω × v) (5)
onde m é a massa do veı́culo e v é a velocidade do seu baricentro, medida no
Referencial Inercial, mas decomposta sobre os eixos do Referencial Móvel. v é
traduzida como uma combinação linear dos versores de F ′, justifica a aceleração
medida no Inercial, v̇, conter o termo ω× v, já que essa base de vetores unitários
muda de direção com o tempo através da rotação ω:
v =
3∑
i=1
v′iê
′
i (6)
F =
3∑
i=1
F ′i ê
′
i (7)
Já a dinâmica de rotação, representada por 4, em função da velocidade de
rotação do veı́culo, assume a seguinte forma:
N = J · ω̇+ ω× (J · ω) (8)
onde ω é a velocidade angular do avião em torno do centro de massa e J é o
tensor de inércia baricêntrico do avião. Assim como no caso da translação, essas
quantidades estão representadas por suas componentes no sistema móvel:
ω =
3∑
i=1
ω′iê
′
i (9)
J =
3∑
i=1
3∑
j=1
J ′ijê
′
i⊗ ê′j (10)
N =
3∑
i=1
N ′i ê
′
i (11)
Adicionalmente, a formulação do problema requer uma representação para
os graus de liberdade da aeronave que possibilite uma descrição unı́voca do seu
movimento. Por isso, foram escolhidos os Ângulos de Euler na sequência 321 para
6
a rotação e as coordenadas cartesianas do Referencial Inercial para a translação,
coordenadas essas conectadas às equações da Dinâmica através das seguintes
relações cinemáticas:
Θ̇ = M(Θ) ω (12)
F ẋ = [R(Θ)]−1 F
′
ẋ = [R(Θ)]−1 v (13)
Θ é a matriz coluna que contém os Ângulos de Euler e R(Θ) a matriz de
rotação que transforma as componentes de vetores no Referencial Inercial nas
suas componentes sobre o Referencial do Corpo. Ainda, F ẋ é a velocidade da
aeronave decomposta no Sistema Inercial, enquanto F ′ẋ é o mesmo vetor veloci-
dade, porém expresso no Sistema do Corpo e já definido anteriormente como v.
M(Θ) é uma matriz cujas entradas são funções dos ângulos de Euler dadas na
sequência com as demais definições aqui mencionadas:
F ẋ =
3∑
i=1
ẋiêi (14)
Θ =
φθ
ψ
 (15)
R(Θ) =
1 0 00 cosφ − senφ
0 senφ cosφ
 cos θ 0 senθ0 1 0
− senθ 0 cos θ
cosψ − senψ 0senψ cosψ 0
0 0 1

(16)
M(Θ) =
1 0 − senθ0 cosφ cos θ senφ
0 − senθ cos θ cosψ
 (17)
De forma suscinta, o modelo matemático que possibilita, juntamente com as
condições iniciais, fazer previsões acerca da evolução temporal dos graus de liber-
dade da aeronave é composto por 5, 8, 13 e 12. Porém, cabe rearranjá-las para
isolar as derivadas de maior ordem em um vetor, enquanto o restante dos termos
é agrupado em outro vetor da seguinte maneira:
η̇ = f(η,U , t) (18)
7
onde η é a matriz coluna que contém os graus de liberdade e as velocidades a
eles associadas, U é o vetor que comporta as entradas no sistema (neste caso, as
resultantes da força e torque externos) e f(η,U , t) é a matriz coluna que carrega
as não-linearidades das equações diferenciais:
η =

Fx
Θ
v̇v̂
ω
 (19)
U =
(
F
N
)
(20)
f(η,U , t) =

[R(Θ)]−1 v
M(Θ) ω
m−1F − ω× v
J−1[N − ω× (J · ω)]
 (21)
2.3 Dinâmica de Voo
Ao aplicar as equações da Dinâmica do Corpo Rı́gido ao avião, os esforços a
serem levados em conta são as forças e torques aerodinâmicos, a gravidade e as
forças e torques propulsivos (cuja propagação se inicia com as turbinas).
O modelo dos esforços aerodinâmicos depende, fundamentalmente, do regime
de escoamento do fluido no qual a aeronave está imersa e das propriedades geométricas
do veı́culo, uma vez que as interações fluido-estrutura se dão por forças de su-
perfı́cie. Tais caracterı́sticas aerodinâmicas se traduzem em um rol de variáveis
que vão desde parâmetros da cinemática do sistema aeronave-fluido às deflexões
de regiões móveis do veı́culo, sendo estas últimas controláveis pelo piloto au-
tomático do avião.
F = Fg + Fp + Fa (22)
N = Np +Na (23)
onde Fg é a força peso, Fp é a força de propulsão da hélice, Fa é a força
aerodinâmica,Na é o torque aerodinâmico eNp o torque propulsivo, dados pelas
8
equações que seguem:
Fg = R(Θ)
 00
mg
 (24)
Fp =
1
2
ρSpropCprop
(kmotoruo)2 − V 2a0
0
 (25)
Fa =
1
2
ρSV 2a
 C1(α) + Cω21 c2Vaω2 + Cue1 (α)ueC02 + Cβ2 β + (Cω12 ω1 + Cω32 ω3) b2Va + Cua2 (α)ua + Cur2 (α)ur
C3(α) + C
ω2
3
c
2Va
ω2 + C
ue
3 (α)ue

(26)
Np =
1
2
ρSpropCprop
−kT (kPuo)20
0
 (27)
Na =
1
2
ρSV 2a

b(C0l + C
β
l β + (C
ω1
l ω1 + C
ω3
l ω3)
b
2Va
+ Cual (α)ua + C
ur
l (α)ur)
c(C0m + C
α
mα + C
ω2
m
c
2Va
ω2 + C
ue
m ue)
b(C0n + C
β
nβ + (C
ω1
n ω1 + C
ω3
n ω3)
b
2Va
+ Cuan (α)ua + C
ur
n (α)ur)

(28)
onde os coeficientes dependentes de α são dados por:
C1(α)
Cω21 (α)
Cuel (α)
C3(α)
Cω23 (α)
Cue3 (α)
 =

−(aD + bDα) (aL + bLα)
−Cω2D C
ω2
L
−CueD C
ue
L
−(aL + bLα) −(aD + bDα)
−Cω2L C
ω2
D
−CueL C
ue
D

(
cosα
senα
)
(29)
À excessão de Va, α e β, que variam no tempo, os demais termos são cons-
tantes que dependem de propriedades geométricas e fı́sicas do ar, do avião e da
interação aerodinâmica destes sistemas, assumindo escoamento laminar desen-
volvido e aproximações lineares de coeficientes de arrasto e sustentação com α.
Uma discussão detalhada deste tópico pode ser encontrada em [1].
9
Figura 3: Velocidade da aeronave em relação ao ar.
Va especifica a norma do vetor velocidade da aeronave medida no suporte
material por onde ela trafega, enquanto α e β dão a orientação deste vetor em
relação ao Referencial do Corpo, representado na Figura 3 por Va.
Incorporando os modelos das forças e torques mostrados nesta sessão no vetor
U , presente na Equação 20 da sessão anterior, e considerando que Va também
pode ser especificado em função dos graus de liberdade da aeronave, chega-se ao
mesmo conjunto de equações diferenciais, porém com a substituição do vetor de
forças:
η̇ = f(η,u, t) (30)
u, neste caso, é o vetor de deflexões das superfı́cies de controle, dado por:
u =

ue
uo
ua
ur
 (31)
10
ue é a deflexão do profundor, uo é um número adimensional que representa a
porcentual da máxima capacidade de rotação do eixo do motor da hélice requisita-
do nas manobras, ua é a deflexão do aileron e ur a orientação do leme da aeronave.
Tais componentes do vetor u, juntamente com os atuadores a eles associados são
ilustrados na figura abaixo:
Figura 4: Atuadores na aeronave e esforço de controle a elas associado - diagrama
adaptado de [1].
Vale ressaltar que uo é proporcional à rotação Ω do motor, conforme mostra a
Figura 4.
2.4 Trimagem e Linearização do Modelo
Antes de iniciar a formulação das leis de controle, é necessário achar uma
região de operação da aeronave, ou seja, uma condição de vôo em torno da qual ela
consiga permanecer estável. Isso implica em determinar a forma como o veı́culo
deve se movimentar e calcular quais os comandos devem ser inferidos ao mesmo
para a execução do vôo escolhido.
Ou seja, a trimagem envolve a escolha de pontos η∗ e η̇∗ que se traduzam no
vôo de operação e na busca de um vetor de comando u∗ que produza os estados
e as derivadas temporais dos estados desejados. Então, o problema se reduz à
resolução do seguinte sistema de equações não-lineares algébricas:
11
η̇∗ = f(η∗,u∗, t) (32)
onde u∗ é a solução do sistema de equações 32 para uma dada trajetória(η∗, η̇∗), em regime permanente.
Essa solução não é simples de ser encontrada, uma vez que não há resposta
analı́tica para muitos sistemas não-lineares, incluindo o de aeronaves rı́gidas. E a
resolução do problema envolve técnicas de otimização, em particular, algoritmos
de programação não linear. Dentro dessa classe de algoritmos, foi utilizada a
Programação Sequencial Quadrática, implementada em uma função MATLAB
conhecida por trim function, que é chamada em um código que será mostrado no
fim desta sessão.
Além do conhecimento do ponto de trimagem, sobre o qual o vôo ocorrerá, é
necessário fazer previsões acerca da resposta temporal do sistema para situações
onde novas forças e torques possam vir a solicitá-lo. Porém, além da própria
resolução de um sistema algébrico não-linear ser onerosa e inviável, conforme já
discutido, a integração de equações diferenciais ordinárias não-lineares, do tipo
ẋ = f(x,u, t), guarda dificuldades de cálculo ainda maiores de um ponto de
vista analı́tico. Por isso, é interessante tentar estratégias alternativas, inclusive,
à otimização não linear, já que esta última demandaria esforço computacional
mais alto no cálculo de respostas transitórias do que no de resposta em regime
permanente.
Uma saı́da pertinente é a utilização da Teoria da Perturbação, que pressupõe
assumir um ponto de operação, Γ∗, previamente calculado na trimagem, e admite
que a evolução temporal do sistema ocorre nas suas vizinhanças, ou seja, os graus
de liberdade variam em pequenas perturbações contabilizadas a partir do ponto
de operação. O range dessa variação é conhecido como região de operação e
denotado porRδ:
Γ∗ = (u∗,η∗, η̇∗) (33)
Rδ = (u∗ + δu,η∗ + δη, η̇∗ + δη̇) (34)
Avaliando, então, as equações de movimento 30 na região de operação, Rδ,
podemos reescrevê-la através de diferenciais inexatas, ou seja, com o operador
δ(·) aplicado à esquerda e à direita da equação:
δη̇ = δf(η∗ + δη,u∗ + δu, t) (35)
12
δf , à esquerda de 35, pode ser expresso da seguinte forma:
δη̇ = A δη +B δu (36)
onde
A =
∂f
∂η
∣∣∣∣
Γ∗
(37)
B =
∂f
∂u
∣∣∣∣
Γ∗
(38)
Para cada regime de voo da aeronave, além do par (A,B), é conveniente,
muitas vezes, definir um vetor ν que contenha as variáveis que se quer controlar
durante a operação do veı́culo.
Normalmente, as componentes de δν são combinações lineares das compo-
nentes de δη e esses vetores se conectam por uma matriz, denotada porC, definida
pela seleção das variáveis de controle envolvidas no projeto.
δν = C δη (39)
A utilização de modelos lineares descritos pelas equações 36 e 39 será eluci-
dada na sessão que segue, que contará com a escolha de uma região de operação
e com o cálculo de modelos de função de transferência para a aeronave.
2.5 Funções de Transferência de Regimes de Voo: O Caso da
Curva Coordenada
Para aeronaves comerciais e mesmo para alguns veı́culos militares, é possı́vel
assumir a dinâmica de voo em dois modos distintos: o modo latero-direcional e o
longitudinal.
O modo latero-direcional diz respeito aos graus de liberdade associados à ro-
lagem, à guinada e à translação, sendo esta última no eixo x′1, enquanto o longi-
tudinal se caracteriza pela arfagem e por duas translações, uma em x′2 e outra em
x′3.
A dinâmica linear desses modos será aqui explorada em um regime de voo
conhecido por Curva Coordenada. Nela, a aeronave circunda um alvo, do qual
mantém uma distância R, aumenta sua altitude, fazendo o segmento de reta do
seu maior comprimento, da pompa à proa, se inclinar de um ângulo γ a partir do
plano x1x2. R e γ são mostrados na imagem que segue:
13
Figura 5: Aeronave e alvo na Curva Coordenada.
Sob a hipótese de que o vento se encontra estacionário em relação ao Refe-
rencial Inercial, ou seja, v = Va, assumindo ẋ3 constante e que a aceleração do
centro de massa aponta na direção x′1 × x3, é possı́vel definir um vetor η̇∗ da
seguinte forma:
η̇∗ =

valor qualquer
outro valor qualquer
V ∗a senγ
∗
05×1
V ∗a
R∗
cos γ∗
03×1

(40)
R∗, V ∗a e γ
∗ caracterizam de forma unı́voca a trimagem. E através da Equação 32,
com o valor de η̇∗, são encontrados os valores de η∗ e u∗, que constituem o ponto
de operação Γ∗, ponto este determinante para o cálculo das matrizes A e B por
37 e 38, respectivamente.
No entanto, as variáveis de interesse ao controle dessa manobra compõe um
14
subconjunto ν das variáveis de estado, cuja seleção se dá pela Equação 39. C e
δν são enunciados na sequência:
C =

01×2 −1 01×9
01×3 1 01×8
01×7
1
V ∗a
01×4
01×3
Rg(secφ∗)2√
tgφ∗
01×8
 (41)
δν =

δh
δφ
δβ
δVa
 (42)
h é a altitude atingida pela aeronave, dada pelo negativo da coordenada 3 do
Referencial Inercial:
h = −x3 (43)
Um modelo do sistema em funções de transferência, que expressa direta-
mente as relações saı́da/entrada, pode ser extraı́do a partir das Equações 36 e
39, operando ambas as igualdades com a Transformada de Laplace, L(·):
G = C(sI −A)−1 (44)
Chega-se, portanto, a um modelo pertinente ao projeto de compensadores
clássicos, cujas técnicas serão exploradas na sessão que segue.
Código matlab 1: Cálculo do modelo linear - programa adaptado de [2].
1 P . g r a v i t y = 9 . 8 ;
2
3 %p h y s i c a l p a r a m e t e r s o f a i r f r a m e
4 P . mass = 1 . 5 6 ;
5 P . Jx = 0 . 1 1 4 7 ;
6 P . Jy = 0 . 0 5 7 6 ;
7 P . Jz = 0 . 1 7 1 2 ;
8 P . Jxz = 0 . 0 0 1 5 ;
9
10 % aerodynamic c o e f f i c i e n t s
11 P .M = 5 0 ;
15
12 P . e p s i l o n = 0 . 1 5 9 2 ;
13 P . a l p h a 0 = 0 . 4 7 1 2 ;
14 P . rho = 1 . 2 6 8 2 ;
15 P . c = 0 . 3 3 0 2 ;
16 P . b = 1 . 4 2 2 4 ;
17 P . S wing = 0 . 2 5 8 9 ;
18 P . S prop = 0 . 0 3 1 4 ;
19 P . k moto r = 2 0 ;
20 P . C L 0 = 0 . 2 8 ;
21 P . C L a lpha = 3 . 4 5 ;
22 P . C L q = 0 . 0 ;
23 P . C L d e l t a e = −0.36;
24 P . C D 0 = 0 . 0 3 ;
25 P . C D q = 0 . 0 ;
26 P . C D d e l t a e = 0 . 0 ;
27 P . C M 0 = 0 . 0 ;
28 P . C M alpha = −0.38;
29 P . C M q = −3.6;
30 P . C M d e l t a e = −0.5;
31 P . C Y 0 = 0 . 0 ;
32 P . C Y beta = −0.98;
33 P . C Y p = −0.26;
34 P . C Y r = 0 . 0 ;
35 P . C Y d e l t a a = 0 . 0 ;
36 P . C Y d e l t a r = −0.17;
37 P . C e l l 0 = 0 . 0 ;
38 P . C e l l b e t a = −0.12;
39 P . C e l l p = −0.26;
40 P . C e l l r = 0 . 1 4 ;
41 P . C e l l d e l t a a = 0 . 0 8 ;
42 P . C e l l d e l t a r = 0 . 1 0 5 ;
43 P . C n 0 = 0 . 0 ;
44 P . C n b e t a = 0 . 2 5 ;
45 P . C n p = 0 . 0 2 2 ;
46 P . C n r = −0.35;
47 P . C n d e l t a a = 0 . 0 6 ;
48 P . C n d e l t a r = −0.032;
49 P . C prop = 1 ;
16
50
51
52 % wind p a r a m e t e r s
53 P . wind n = 0 ;
54 P . wind e = 0 ;
55 P . wind d = 0 ;
56 P . L wx = 1250 ;
57 P . L wy = 1750 ;
58 P . L wz = 1750 ;
59 P . sigma wx = 1 ;
60 P . sigma wy = 1 ;
61 P . sigma wz = 0 . 7 ;
62 % P . sigma wx = 0;
63 % P . sigma wy = 0;
64 % P . s igma wz = 0;
65 P . Va0 = 1 3 ;
66
67 % a u t o p i l o t sample r a t e
68 P . Ts = 0 . 0 1 ;
69
70 % compute t r i m c o n d i t i o n s u s i n g ’ m a v s i m c h a p 5 t r i m . mdl
’
71 P . Va = 1 7 ;
72 gamma = 5∗ pi / 1 8 0 ; % d e s i r e d f l i g h t pa th a n g l e (
r a d i a n s )
73 R = 150 ; % d e s i r e d r a d i u s (m) − use (+)
f o r r i g h t handed o r b i t ,
74 % (−) f o r
l e f t handed o r b i t
75 % f i r s t c u t a t i n i t i a l c o n d i t i o n s
76 P . pn0 = 0 ; % i n i t i a l Nor th p o s i t i o n
77 P . pe0 = 0 ; % i n i t i a l Eas t p o s i t i o n
78 P . pd0 = 0 ; % i n i t i a l Down p o s i t i o n ( n e g a t i v e
a l t i t u d e )
79 P . u0 = P . Va ; % i n i t i a l v e l o c i t y a long body x−a x i s
80 P . v0 = 0 ; % i n i t i a l v e l o c i t y a long body y−a x i s
81 P . w0 = 0 ; % i n i t i a l v e l o c i t y a long body z−a x i s
82 P . ph i0 = 0 ; % i n i t i a l r o l l a n g l e
17
83 P . t h e t a 0 = 0 ; % i n i t i a l p i t c h a n g l e
84 P . p s i 0 = 0 ; %i n i t i a l head ing a n g l e
85 P . p0 = 0 ; % i n i t i a l body frame r o l l r a t e
86 P . q0 = 0 ; % i n i t i a l body frame p i t c h r a t e
87 P . r0 = 0 ; % i n i t i a l body frame yaw r a t e
88
89 % run t r i m commands
90 [ x t r i m , u t r i m ]= c o m p u t e t r i m ( ’ mavs im t r im ’ , P . Va , gamma
, R) ;
91 P . u t r i m = u t r i m ;
92 P . x t r i m = x t r i m ;
93
94 % s e t i n i t i a l c o n d i t i o n s t o t r i m c o n d i t i o n s
95 % i n i t i a l c o n d i t i o n s
96 P . u0 = x t r i m ( 4 ) ; % i n i t i a l v e l o c i t y a long body x
−a x i s
97 P . v0 = x t r i m ( 5 ) ; % i n i t i a l v e l o c i t y a long body y
−a x i s
98 P . w0 = x t r i m ( 6 ) ; % i n i t i a l v e l o c i t y a long body z
−a x i s
99 P . ph i0 = x t r i m ( 7 ) ; % i n i t i a l r o l l a n g l e
100 P . t h e t a 0 = x t r i m ( 8 ) ; % i n i t i a l p i t c h a n g l e
101 P . p0 = x t r i m ( 1 0 ) ; % i n i t i a l body frame r o l l r a t e
102 P . q0 = x t r i m ( 1 1 ) ; % i n i t i a l body frame p i t c h
r a t e
103 P . r0 = x t r i m ( 1 2 ) ; % i n i t i a l body frame yaw r a t e
104
105 % compute d i f f e r e n t t r a n s f e r f u n c t i o n s
106 [ T p h i d e l t a a , T c h i p h i , T t h e t a d e l t a e , T h t h e t a ,
T h Va , T V a d e l t a t , T V a t h e t a , T v d e l t a r ] . . .
107 = c o m p u t e t f m o d e l ( x t r i m , u t r i m , P ) ;
108
109 %l a t
110 [ n ph i , d p h i ]= t f d a t a ( T p h i d e l t a a , ’ v ’ ) ;
111 [ n v d e l r , d v d e l r ]= t f d a t a ( T v d e l t a r , ’ v ’ ) ;
112
113 %l o n
114 [ n t h e t a , d t h e t a ]= t f d a t a ( T t h e t a d e l t a e , ’ v ’ ) ;
18
115 [ n h t h , d h t h ]= t f d a t a ( T h t h e t a , ’ v ’ ) ;
116 [ n Va th , d V a t h ]= t f d a t a ( T V a t h e t a , ’ v ’ ) ;
117 [ n V a d e l t , d V a d e l t ]= t f d a t a ( T V a d e l t a t , ’ v ’ ) ;
118
119
120 % l i n e a r i z e t h e e q u a t i o n s o f mot ion around t r i m
c o n d i t i o n s
121 [ A lon , B lon , A l a t , B l a t ] = c o m p u t e s s m o d e l ( ’
mavs im t r im ’ , x t r i m , u t r i m ) ;
3 Projeto de Compensadores de Voo
O controle de voo, definido para a curva coordenada, é projetado para os mo-
dos latero-direcional e longitudinal separadamente. E as funções de transferência
que constam emG valem para a seguinte terna de parâmetros:
(V ∗a , γ
∗, R∗) = (17m/s,
5π
180
rad, 150m) (45)
A metodologia de projeto dependerá da estrutura algébrica das funções de
transferência da diagonal de G e dos acoplamentos existentes entre as variáveis
que formam o vetor δν, acoplamentos estes dados pelos elementos fora dessa dia-
gonal. As plantas Gij , com os valores numéricos dos coeficientes dos polinômios
que as compõe, serão mostradas a seguir, juntamente com o procedimentos de
sı́ntese dos compensadores.
O controle será feito em malha fechada, estratégia justificada pelo fato de que
aeronaves devem ser robustas a distúrbios, sejam eles ambientais, como vento, ou
paramétricos, como variações de massa do veı́culo, por exemplo.
3.1 Controle do voo latero-direcional
No modo latero-direcional, as únicas variáveis a serem controladas são a ro-
lagem, φ(t), e o ângulo de derrapagem, β(t).
Como essas plantas são desacopladas, serão explorados conceitos adicionais
de robustez a perturbações na resolução de seus problemas.
19
3.1.1 Rolagem
Na rolagem, o esforço de controle que a gera diretamente é a deflexão do
aileron, δua(t), e a função de transferência que conecta essas duas grandezas, da
perspectiva de entrada (aileron) e saı́da (rolagem), é a função G22 da matriz G,
dada, a partir da trimagem, pela seguinte razão polinomial:
δφ(s)
δua(s)
= G22(s) =
25,83
s(s+ 4,72)
(46)
Deseja-se encontrar um compensador C2(s) que infira deflexões no aileron
para que a planta responda com uma rolagem comandada pelo piloto, φc(t). Ou
seja, C2(s) deve garantir a validade da igualdade δφ(t) = δφc(t).
Figura 6: Sistema de controle de rolagem.
Porém, nota-se que G22(s) é uma razão polinomial da seguinte forma:
G(s) =
ω2n
s(s+ 2ζωn)
(47)
Figura 7: Sistema dinâmico linear de segunda ordem.
20
A frequência caracterı́stica ou natural, ωn, é uma métrica da rapidez com que
as oscilações do sistema ocorrem na ausência de esforços aplicados, enquanto o
fator de amortecimento, ζ , exprime a taxa de dissipação da energia cinética para
o ambiente, ou seja, maneira pela qual as vibrações são atenuadas.
Uma função de transferência que possa ser representada através de um sistema
em malha fechada com G(s) na malha aberta, conforme mostrado na 7, é um
sistema de segunda ordem que exibe uma resposta similar a do gráfico abaixo
quando é comandada por um sinal yc = 1:
Figura 8: Resposta de um sistema de segunda ordem ao degrau unitário.
A resposta de
G
1 +G
converge segundo o padrão da Figura 8 por duas razões:
• G possui um polo na origem, o que possibilita erro nulo ao degrau unitário
em regime estacionário, y(∞)− yc = 0 ;
• os valores singulares da razão y(s)
yc(s)
possuem parte real negativa, conforme
mostra a Figura 9, e implicam em uma resposta y(t) modulada por uma
função exponencial de argumento negativo, e−ζωnt, ou seja, resultam em
uma resposta limitada;
A função de transferência do sistema em malha fechada, descrito pela Figura
7, é dada pela seguinte equação:
21
y(s)
yc(s)
=
G(s)
1 +G(s)
=
ω2n
s2 + 2ζωns+ ω2n
(48)
A equação caracterı́stica desse sistema acima é dada pela nulidade do deno-
minador, ou seja, é uma equação que expõe as singularidades da função
y(s)
yc(s)
:
1 +G(s∗) = 0 (49)
Os valores singulares de
y(s)
yc(s)
(no caso, um sistema de segunda ordem),
também chamados de pólos do sistema, são as raı́zes da equação 49, ou seja,
os valores de s∗ que satisfazem 1 +G = 0:
s∗ = −ζωn ± iωn
√
1− ζ2 = ωne±i[π−arccos(ζ)] (50)
O lugar geométrico dessas raı́zes no Plano de Argand-Gauss é mostrado a
seguir:
Figura 9: Lugar das raı́zes de um sistema de segunda ordem amortecido.
22
No caso, os pólos se encontram no semiplano esquerdo do Diagrama de Ar-
gand que, por se tratar de um plano em que são representados valores da variável
s, é também conhecido, como plano s.
A localização dos pólos no plano s determina se a resposta do sistema a uma
entrada teste convergirá ou se será desviada de forma não limitada.
Figura 10: Análise de estabilidade do sistema em malha fechada.
O eixo imaginário do Diagrama de Argand é a fronteira que separa a região
de estabilidade da zona de instabilidade. E para o caso de pólos sobre este eixo,
Re(s∗) = 0, o sistema responde a uma entrada degrau, yc = 1, com uma oscilação
não-amortecida em torno da entrada. Entretanto, a saı́da não converge para ne-
nhum valor particular, tampouco diverge:
23
Figura 11: Eixo imaginário como limite de estabilidade do sistema.
Em regiões fora do eixo imaginário, os pólos implicam ou em respostas limi-
tadas e convergentes, para Re(s∗) < 0, ou em respostas que crescem de maneira
indiscriminada no tempo, no caso, quandoRe(s∗) > 0.
24
Figura 12: Conexão entre convergência nas respostas e pólos de G
1+G
.
E como G22(s) possui similaridade com o formato padrão de um sistema de
segunda ordem, representado por G, será utilizado um compensador C2 estático e
de valor unitário:
C2(s) = 1 (51)
Assim, a função de transferência de malha aberta, oriunda da equivalência
entre C2G22 e G nas figuras 6 e 7, assume a seguinte forma:
C2G22(s) =
25,83
s(s+ 4,72)
=
ω2n
s(s+ 2ζωn)
(52)
Por inspeção, com base em 52, é previsto que o sistema C2G22
1+C2G22
seguirá sinais
do tipo degrau (funções de valor constante no tempo) com erro nulo (pelo pólo
na origem de C2G22) e com amortecimento e frequência natural assumindo os
valores presentes no par (ωn; ζ) = (5,08rad/s; 0,46).
25
A resposta do sistema δφ(s)
δφc(s)
, φ(t), a uma entrada degrau de 1◦, φc(t)= π180rad,
foi plotada com a ajuda do MATLAB:
Figura 13: Resposta de δφ(s)
δφc(s)
a uma entrada degrau de π
180
rad.
Apesar da resposta satisfatória, o sistema precisa ser capaz de rejeitar perturbações
não modeladas. No diagrama de blocos que segue, é exibida a ação de um distúrbio
na planta:
Figura 14: Sistema de controle de δφ com um distúrbio na planta.
δφ(s) =
[
C2G22(s)
1 + C2G22(s)
]
δφc(s) +
[
G22(s)
1 + C2G22(s)
]
d(s) (53)
Na Equação 53, verifica-se que o efeito do distúrbio d(s) é minimizado e que
26
a saı́da δφ(s) rastreia o comando δφc(s) se o compensador assumir valores altos
segundo a comparação contida na proposição abaixo:
[C2G22(s) ≫ 1]→ δφ ≈ δφc (54)
Uma simulação de um distúrbio foi executada no Simulink. Um sinal degrau
somado a 4 funções harmônicas foi tomado como distúrbio conforme o diagrama
abaixo:
Figura 15: Simulação simulink da atuação de C2 = 1
frente ao distúrbio d(t)
d(t) =
π
180
[5 + 4 sen(t) + 3 sen(2t) + 2 sen(3t) + sen(4t)] rad (55)
Mantendo o mesmo compensador, a resposta do sistema na presença do distúrbio
mostra que sua capacidade de seguir o comando foi deteriorada:
27
Figura 16: Gráfico da resposta y(t) sob efeito de um distúrbio na planta.
E para diminuir a importância do distúrbio sobre y, aumentaremos o valor de
C2(s), mantendo o mesmo ainda como uma constante. A previsão, por análise da
proposição 54, é de que a resposta se aproxime mais do comando yc:
C2 = 100 (56)
28
Figura 17: Efeito do aumento do ganho C2 na resposta.
Na Figura 17, verifica-se que o efeito do distúrbio diminuiu. Porém, o aumento
do ganho do compensador implicou em oscilações de alta frequência da resposta
de regime transiente e em valores que excedem o dobro do valor do comando dado
ao sistema.
Uma maneira de contornar este problema é manter o valor do compensador
da tentativa inicial, C2 = 1, e construir um sub-sistema auxiliar que seja capaz
de estimar o efeito do distúrbio e calcular um esforço de controle adicional que
se contraponha ao mesmo. Um arranjo possı́vel para a execução dessa lógica de
controle é dada no diagrama de blocos abaixo:
29
Figura 18: Estimador e regulador de distúrbio incorporados ao sistema de contro-
le.
A estimação do efeito do distúrbio consiste em subtrair o efeito de δua sobre a
resposta do seu efeito combinado ao do distúrbio, resultando em uma observação
isolada do efeito do distúrbio sobre a saı́da.
O regulador, por sua vez, compara o efeito do distúrbio com o valor 0, como
um comando de minimização do efeito do distúrbio, e multiplica o erro por um
ganho K, resultando em um sinal de controle que tenderá à função −d(t) para o
cancelamento da referida perturbação.
30
Figura 19: Compensador, estimador e regulador em ambiente Simulink.
Foi escolhido um valor alto para o ganho do regulador para que o efeito do
distúrbio fosse diminuı́do significativamente. A simulação, em ambiente Simulink,
resultou na resposta mostrada no gráfico que segue:
K = 10000 (57)
Figura 20: Resposta do sistema depois da inclusão do regulador/estimador.
31
Com isso, finaliza-se o projeto do compensador de rolagem. A estrutura pro-
posta, por ter utilizado o modelo da planta na lógica de controle, é conhecida
como controlador com modelo interno. Esta classe de controladores não será ex-
plorada em maiores detalhes, já que o objetivo aqui é explorar estritamente con-
ceitos básicos da Teoria de Controle. E por mais que este tipo de estratégia faça
parte de cursos mais avançados de controle, o raciocı́nio em cima da construção
de sua estrutura algébrica ainda é simples. Além disso, o controlador com modelo
interno foi selecionado para a rejeição de distúrbios pelo fato de que métodos al-
ternativos exigiriam do leitor domı́nio pleno de análise da resposta em frequência
da planta, tema que poderá ser melhor apresentado depois que o conceito de Lugar
das Raı́zes for introduzido (ainda na subsessão seguinte).
Maiores detalhes sobre essa abordagem podem ser encontrados em [3].
3.1.2 Ângulo de Derrapagem
De forma similar ao caso da rolagem, o sistema de controle do ângulo de
derrapagem pode ser descrito por um diagrama de blocos segundo a figura abaixo:
Figura 21: Sistema de controle do ângulo de derrapagem.
δβc é o comando, δβ o ângulo de derrapagem, C3 o compensador a ser proje-
tado, δur a deflexão do leme de direção e G33 a função de transferência da planta,
dada a seguir:
δβ(s)
δur(s)
= G33(s) =
−5,17
s+ 1,75
(58)
A Transfordada de Fourier de G33(s), G33(s = wi), revela que o ganho que o
sistema confere às entradas de frequências muito baixas (w ≈ 0) é negativo. Ou
seja, para sinais de frequência nula como a função degrau, por exemplo, a resposta
do sistema em regime permanente terá o sinal oposto ao da entrada.
32
G33(0) = −
5,17
1,75
(59)
Como é de interesse que o sistema rastreie funções e preserve o sinal dos
valores que ela assume no tempo, o compensador C33 terá de incorporar o valor
−1 entre os fatores que vierem a compor sua estrutura.
Além disso, outra exigência de desempenho tı́pica que a planta deve atender,
neste caso, é de exibir capacidade de seguir uma função de derivada não nula. Um
exemplo de comando tı́pico é para que a saı́da atinja um determinado valor por
uma função que cresce segundo uma reta até o instante que o alcança:
βc(t) =

0,25t se t ≤ 4
1 se t > 4
(60)
Como parte do sinal é uma função rampa, definida no intervalo [0; 10s], a
função de transferência de malha aberta precisa de dois pólos na origem para
rastrear o comando com um erro nulo. Já que G33 não possui essa propriedade,
C3 deve comportá-los em sua estrutura. Por isso, C3 deve ser definido segundo
uma razão de polinômios que contenha−1 e 1/s2 como fatores que o constituem:
C3(s) = −
k
s2
(61)
No entanto, a adição de pólos na malha aberta implica, normalmente, no deslo-
camento dos pólos de malha fechada para o semi-plano direito do Diagrama de
Argand, ou seja, para a região de instabilidade. E para verificar a localização
destes pólos de malha fechada, basta variar o ganho k de 0 a ∞ e achar uma
solução da equação caracterı́stica do sistema para cada compensador encontrado
sua excursão:
1 + C3G33(sk) = 0 (62)
sk são os pólos da função de transferência de malha fechada para um dado
valor de k. Reescrevendo a equação caracterı́stica para o compensador e a planta
especificados, tem-se que:
1 +
k
s2k
5,17
(sk + 1,75)
= 0 (63)
Para cada valor de k, neste caso, existem 3 pólos a serem marcados no Di-
agrama de Argand. Em k = 0, os pólos de malha fechada assumem posições
33
que coincidem com os pólos de malha aberta: sk=0 = 0; 0;−1,75. E se o ganho
variar somente por valores positivos, os pólos de malha fechada assumem outras
posições. Os afixos correspondentes a estas últimas são encontrados através da
Equação 63 na iteração de k no intervalo [0;∞).
Figura 22: Variação do posicionamento dos polos em malha fechada à medida que
o ganho k é excursionado.
O Diagrama de Argand é comumente chamado de Lugar das Raı́zes quando
é utilizado para mostrar a conexão entre as diversas posições assumidas pelos
pólos de uma função de transferência com os possı́veis valores de uma determi-
nada propriedade daquele sistema (pode ser uma massa, uma resistência ou, como
mostrado até agora, um ganho do controlador).
Na Figura 22, observa-se que quaisquer valores de k levam 2 dos 3 pólos de
malha fechada a se posicionarem no semiplano direito do plano s, evidenciando
que C3 leva o sistema em malha fechada à instabilidade. Por isso, a estrutura
do compensador deve ser modificada, visando um percurso estável destes pólos
atrelado à variação dos parâmetros livres do controlador.
Uma técnica convencional de mudar a localização dos pólos de malha fechada
é inserir zeros na função de transferência de malha aberta. Zeros são os valoresde s que zeram o numerador de uma função de transferência; para H(s) = (s +
1)/(s + 10), por exemplo, o zero é o valor s = −1, já que ele faz com que H(s)
se anule.
Os zeros de malha aberta são pontos atratores no Lugar das Raı́zes, ou seja,
à medida que o ganho k da malha aberta aumenta, os pólos de malha fechada
34
modificam suas posições para regiões cada vez mais próximas das coordenadas
desses zeros. E quando o ganho atinge um valor muito alto, k −→ ∞, os pólos
de malha fechada são alocados nas posições dos zeros de malha aberta, porém,
na razão de um pólo na posição de um zero, com os pólos excedentes localizados
em outras regiões. Então, com base nessa atração que os zeros de malha aberta
inferem sobre os pólos de malha fechada, uma tentativa pertinente é incorporar
zeros estáveis em C3.
Inicialmente, será inserido um zero no semiplano esquerdo, cujo efeito será
verificado por meio do Diagrama do Lugar das Raı́zes plotado para a variação de
k. Para isso, C3 será redefinido, multiplicando sua estrutura anterior, exibida na
Equação 61, por um zero estável, escolhido em s = −3 neste caso:
C3(s) = −
k(s+ 3)
s2
(64)
Reescrevendo a equação caracterı́stica, 62, com o segundo C3, presente na
Equação 64, chega-se à seguinte equação caracterı́stica:
1 +
k(sk + 3)
s2k
5,17
(sk + 1,75)
= 0 (65)
Plotando o Diagrama do Lugar das Raı́zes dessa nova equação caracterı́stica,
observa-se o novo caminho descrito pelos pólos de malha fechada na excursão de
k:
35
Figura 23: Novo Diagrama do Lugar das Raı́zes para um novo C3.
Pela figura acima, observa-se que o caminho dos pólos de malha fechada se
modificou pela ação do zero adicionado, porém, a adição desse elemento não foi
o bastante para atrair os ramos que se propagavam para o semiplano direito. Neste
caso, a segunda tentativa pertinente é a adição de outro zero estável.
O segundo zero será escolhido na mesma posição do primeiro (sua localização
pode ser qualquer desde que seja na região estável do locus), levando C3 a ser
redefinido segundo a seguinte equação:
C3(s) = −
k(s+ 3)2
s2
(66)
Enunciando novamente a equação caracterı́stica do sistema, a Equação 64,
mas com esta última atualização do compensador, chega-se à seguinte igualdade:
1 +
k(sk + 3)
2
s2k
5,17
(sk + 1,75)
= 0 (67)
A plotagem do Lugar das Raı́zes desta última equação caracterı́stica mostra
que os pólos caminham para regiões de resposta amortecida, conforme mostra a
36
figura abaixo:
Figura 24: Pólos de malha fechada estáveis depois da inserção do segundo zero,
para k = 30.
Conforme mostra a Figura 24, foram escolhidos pólos na região de fator de
amortecimento unitário, para que fossem eliminadas oscilações que ocorreriam
caso essas singularidades possuı́ssem partes imaginárias não nulas. O valor do
ganho correspondente a tais pólos é exibido na própria figura 24 e na equação que
segue:
k = 30 (68)
Além destes pólos mencionados, a mesma Figura 24 mostra que um dos ramos
do Lugar das Raı́zes (mostrado pela flecha que vai para a esquerda) evidencia
a existência de um terceiro pólo de coordenada (−150; 0) no plano complexo,
segundo a Figura 25.
37
Figura 25: Pólo em (−150; 0).
Porém, considerando somente a rastreabilidade do sistema, ou seja, a capaci-
dade do mesmo seguir o comando δβc, o peso desse pólo não é relevante, uma vez
que o efeito dele sobre a resposta se dá por uma janela de tempo muito pequena.
Isso ocorre porque os modos do sistema são mais rápidos quanto mais distantes
os pólos a eles associados estiverem do eixo imaginário, conforme mostra a figura
26.
38
Figura 26: Comparação de velocidade dos modos de um sistema dinâmico.
Uma regra prática para atestar se o efeito de um pólo mais rápido pode ser des-
prezado, ou seja, se a velocidade do seu decaimento será alta quando comparada
à rapidez dos demais modos, é verificar se sua distância do eixo imaginário é
maior do que cinco vezes a distância dos outros pólos ao mesmo eixo. E como a
coordenada do pólo mostrado na Figura 25 satisfaz esta regra (comparando com
as coordenadas dos pólos da Figura 24), seu efeito pode ser negligenciado para a
rastreabilidade.
O projeto do sistema de controle do ângulo de derrapagem estaria finalizado
nesta última etapa. Porém, é interessante submeter o sistema a uma perturbação
e verificar se ele é capaz se seguir um comando mesmo na presença dela. Mas,
diferente de um distúrbio sobre a planta, será considerado agora um ruı́do no
sensor, conforme mostra o diagrama a seguir:
39
Figura 27: Sistema de controle do ângulo de derrapagem com ruı́do no sensor.
O ruı́do no sensor foi definido como uma senóide de alta frequência, 100Hz
neste caso:
nβ(t) =
0,2π
180
sen(2π100t) rad (69)
Considerando o ruı́do definido na Equação 69 e o comando dado pela Equação
60, simulamos o sistema em ambiente MATLAB para verificar a resposta δβ(t),
mostrada no gráfico abaixo:
Figura 28: Resposta do sistema comandado por δβ(t) e perturbado por nβ(t).
Pela Figura 28, observa-se que o sistema é sensı́vel ao ruı́do. E para diminuir
tal sensibilidade, é conveniente analisar de que maneira cada entrada influencia a
resposta do sistema e se é possı́vel diminuir o efeito do ruido apenas com algum
ajuste nos parâmetros de C3(s):
40
δβ(s) =
[
C3G33(s)
1 + C3G33(s)
]
δβc(s)−
[
C3G33(s)
1 + C33G33(s)
]
nβ(s) (70)
Os valores de C3 ditam a relevância tanto do comando δβc quanto do ruı́do nβ
na resposta δβ. E pela Equação 70, percebe-se que quanto maior for o valor do
compensador, mais importante será o valor do ruı́do para a resposta do sistema.
Entretanto, o mesmo peso, em módulo, é dado ao comando.
Considerando, em primeira análise, tal igualdade entre os módulos desses pe-
sos, deve-se buscar outra propriedade que implique em efeitos distintos à transformação
dessas entradas. E o detalhe ao qual o projetista deve se ater para a análise de
sinais, neste caso, é a posição relativa dos pólos de malha fechada no Lugar das
Raı́zes.
O pólo mais distante do eixo imaginário, na coordenada (−150; 0), segundo a
Figura 25, é pouco importante em relação à capacidade de rastreamento do sistema
de controle (explicação na Figura 26). Entretanto, um modo rápido como este
contribui para a amplificação de sinais de alta frequência, enquanto modos mais
lentos (por exemplo, os demais da Figura 25, mais próximos do eixo imaginário)
amplificam sinais de baixa frequência e, inclusive, em algumas circunstâncias,
atenuam o efeito de sinais que oscilam mais rapidamente.
De posse deste conhecimento, pode-se mudar a posição dos pólos no root
locus novamente com um ajuste de ganho, por exemplo, para que o sinal δβc se
torne mais importante quando comparado a ηβ . Mas para que isso seja possı́vel, é
necessário saber qual desses sinais de entrada possui componentes de frequências
mais altas; para isso, basta plotar o gráfico do módulo da Transformada de Fourier
deles, uma vez que a mudança de domı́nio do tempo para o da frequência permite
a visualização imediata das amplitudes associadas às harmônicas que constituem
quaisquer sinais.
Antes de traçá-lo, se faz necessário tomar as Transformadas de Fourier dos
sinais ηβ(t) e δβc(t), cujas funções temporais são descritas nas Equações 69 e 60,
respectivamente.
nβ(f) =
∫ ∞
0
e−i2πftnβ(t)dt = 0,2
[
i2πf
(i2πf)2 + (2π100)2
]
(71)
δβc(f) =
∫ ∞
0
e−i2πftδβc(t)dt =
e−4(i2πf) − 1
4(i2πf)2
(72)
Mais detalhes acerca dessa transformada serão explorados no problema de
controle de voo longitudinal da Sessão 3.2. Porém, neste momento, apenas o
41
gráfico do módulo das transformadas de Fourier das entradas das Equações 71 e
72 será suficiente para proceder com a análise:
Figura 29: Módulo das transformadas de Fourier das entradas do sistema.
Percebe-se, pela Figura 29, que o sinal ηβ possui uma componente em 100 Hz
(só uma harmônicaconstitui esse sinal), enquanto δβc detém componentes em
frequências muito mais baixas. Então, se trouxermos o pólo localizado em (-
150;0) no locus para uma coordenada mais próxima ao eixo imaginário, a in-
fluência do sinal que oscila mais lentamente, δβc(t), no caso, será maior e a
relevância do ruı́do deve diminuir drasticamente.
k = 30 foi o ganho escolhido anteriormente para o compensador descrito pela
Equação 66. Mas agora, afim de deixar os pólos mais lentos, podemos atribuir à
mesma constante um valor menor:
k =
1
2
(73)
Com esse novo k, basta avaliar a posição dos pólos de malha fechada no Di-
agrama do Lugar das Raı́zes e averiguar, pela resposta δβ(t), se a relevância do
ruı́do de fato diminuiu:
42
Figura 30: Novas coordenadas assumidas pelos pólos de malha fechada após a
diminuição do ganho.
Figura 31: Resposta do sistema em malha fechada aos efeitos simultâneos do
comando e do ruı́do para o novo ganho k.
Com o último ajuste de ganho, está encerrado o projeto do compensador, cuja
estrutura assume a seguinte forma depois da atualização do ganho k:
43
C3(s) = −
1
2
(s+ 3)2
s2
(74)
3.2 Controle do voo longitudinal
Para controlar a velocidade do ar medida na aeronave, δVa, e a altitude, δh,
serão utilizadas técnicas de controle no domı́nio da frequência. No entanto, rudi-
mentos de Resposta em Frequência serão aqui discutidos para que se possa dar
ênfase ao problema de controle através de conceitos básicos inerentes a essas es-
tratégias.
3.2.1 Rudimentos de Resposta em Frequência
Essas técnicas possuem duas vantagens relevantes frente ao rootlocus:
• podem ser aplicadas diretamente a modelos de plantas derivados de experi-
mentos de identificação.
• permitem inferir boas previsões acerca da sensibilidade do sistema a sinais
variantes no tempo mais elaborados do que as entradas de teste tı́picas con-
sideradas quando o projeto se dá no domı́nio do tempo.
Como o tópico Identificação de Sistemas Dinâmicos é um curso extenso suportado
por uma densa literatura, este texto se focará somente no segundo aspecto: na
questão da sensibilidade dos sistemas à rapidez de oscilação das entradas.
Uma função temporal arbitrária pode ser representada como uma soma in-
finita de harmônicas. Essa série de senóides e cossenóides é conhecida por Soma
Infinita de Fourier ou Série de Fourier, representada pela figura que segue:
44
Figura 32: Série de Fourier de uma função contı́nua no tempo.
Além disso, cada harmônica possui uma representação no Diagrama de Ar-
gand Gauss, mostrada na figura abaixo:
45
Figura 33: Representação complexa de uma harmônica.
Esta representação complexa é conveniente por permitir a exibição de cada
componente senoidal de um sinal arbitrário através de um ponto no plano com-
plexo. Porém, conforme mostra a figura 33, além de contabilizar a frequência da
harmônica, tal assinatura gráfica tem que levar em conta, também, o tempo t.
É mais pertinente, quando se quer avaliar a velocidade de oscilação das com-
ponentes de um sinal, negligenciar o tempo t e cercar somente a variável frequência,
ω, já que esta última é capaz de caracterizar a harmônica univocamente.
Portanto, vale, neste tipo de problema, lançar mão de uma transformação de
sinal que não exponha o tempo, mas que deixe a frequência explı́cita, preservando
a representação do sinal no plano complexo. Essa transformada é conhecida como
Transformada de Fourier e, curiosamente, pode ser definida a partir da Transfor-
mada de Laplace:
F(·) = L(·)|s=jω (75)
Considerando a Equação 75 e que os sistemas fı́sicos, encarados como plantas
em sistemas de controle, lidam com sinais de diferentes componentes nas suas
46
entradas e saı́das, convém, também, aplicar essa transformada para a obtenção de
um modelo da planta na frequência.
Um sistema dinâmico pode ser entendido como um processo que transforma
entradas em saı́das. Em especial, um sistema dinâmico linear estável de pro-
priedades fı́sicas e geométricas constantes (também chamados de sistemas line-
ares invariantes no tempo - LTI - estáveis) é um processo que, em regime per-
manente, transforma a entrada em uma saı́da proporcional à primeira, provendo
a essa saı́da, também, possı́veis atrasos na emissão das respostas, como mostra a
figura a seguir:
Figura 34: Conceito de resposta em frequência para um sistema LTI.
Entretanto, essa transformação não ocorre de qualquer forma. Em regime
permanente, as caracterı́sticas da saı́da diferem das da entrada por propriedades
da planta que residem no seu modelo no domı́nio da frequência. Esse modelo, por
sua vez, é um número complexo dependente da frequência; através dele, é possı́vel
fazer a previsão de como as amplitudes as entradas se modificarão (multiplicadas
pelo módulo desse número na referida frequência) e como as fases dos sinais serão
alteradas (somando-se às mesmas a fase da planta na frequência da entrada).
Essa representação complexa da planta é conhecida como Resposta em Frequência
do sistema e consiste na Transformada de Fourier do sistema fı́sico, obtida a partir
da Transformada de Laplace da planta segundo a equação que segue:
G(jω) = G(s)|s=jω (76)
G(jω) é um número complexo e possui, para ω ∈ R, uma representação faso-
rial no plano complexo, porém, como um fasor variante na frequência:
47
Figura 35: Diagrama de Nyquist de G(jω).
Essa variação do fasor G na frequência ω possui uma representação gráfica
capaz de exibir o contı́nuo de afixos de G no plano complexo em todo o domı́nio.
Esse diagrama é conhecido como Diagrama de Nyquist, ilustrado na Figura 35.
E através dele, é possı́vel fazer previsões da estabilidade e resposta transitória do
sistema em malha fechada que contém G no seu circuito. Este sistema em malha
fechada é representado no diagrama de blocos abaixo:
Figura 36: G inserida no sistema em malha fechada.
A utilização da Curva de Nyquist para fins de análise se dá através do Critério
da Estabilidade de Nyquist, que pode ser enunciado segundo o seguinte excerto:
48
O número de pólos instáveis do sistema em malha fechada é igual a soma
do número de pólos instáveis da planta com o número de voltas no sen-
tido horário que a curva de Nyquist da planta dá ao redor da coordenada
(−1; 0) do plano complexo.
Matematicamente, o Critério de Nyquist pode ser expresso pela equação abaixo:
nf = np + nv (77)
Nela, np é o número de pólos instáveis da planta, nv é o número de voltas que
a curva de Nyquist dá em torno da coordenada (−1; 0) no plano complexo e nf é
a previsão de quantos pólos instáveis o sistema em malha fechada apresentará.
Para esclarecer esse conceito, a previsão de estabilidade se dará sobre um
exemplo que apresenta a seguinte função de transferência:
G(s) = − 100
(s− 1)(s+ 2)(s+ 5)
(78)
De inı́cio, ao olhar pra planta, constata-se que ela detém um pólo instável,
s = +1. Então, np = 1.
Resta traçar o Diagrama de Nyquist para essa G, substituindo s por jω:
49
Figura 37: Diagrama de Nyquist do primeiro exemplo.
Como a curva enlaça a coordenada (−1; 0) no plano complexo no sentido
horário uma vez, nv = 1. Podemos, portanto, efetuar a soma indicada na Equação
77:
nf = np + nv = 1 + 1 = 2 (79)
Então, a previsão para este caso é de que o sistema em malha fechada,
G
1 +G
,
apresentará 2 pólos instáveis. O leitor pode tirar a prova encontrando as raı́zes da
equação caracterı́stica 1+G(s) = 0 e caso efetue a conta de forma correta, deverá
chegar nas seguintes raı́zes: s = −7,28; s = +0,64± 3,46j.
Além da análise de estabilidade do sistema em malha fechada, é possı́vel ex-
trair métricas do seu comportamento durante o regime transitório, mas somente
em sistemas controlados cuja planta possua pólos estáveis, ou seja, pra casos onde
G não possua pólos no semiplano-direito do plano s, quando np = 0:
nf = nv (80)
50
A Equação77 se reduz à Equação 80 nos casos onde o sistema em malha
aberta é estável; e a estabilidade do sistema operando em malha fechada se res-
tringe à avaliação exclusiva do Diagrama de Nyquist. Para trabalhar com essa
classe de sistemas fı́sicos, o problema será raciocinado em cima de uma G parti-
cular:
G(s) =
100
(s+ 1)(s+ 2)(s+ 5)
(81)
G, neste caso, mostrada na equação acima, possui todos os seus pólos na
região estável do plano s. Logo, se a inserirmos em um sistema em malha fechada,
G
1 +G
, a resposta deste último a um comando será ilimitada somente se a Curva
de Nyquist de G enlaçar a coordenada (−1, 0), como mostra a equação 80.
Traçando, então, o referido Diagrama de Nyquist, chega-se ao seguinte retrato:
Figura 38: Diagrama de Nyquist da G.
Se a mesma curva de Nyquist da figura acima não enlaçasse a coordenada
(−1; 0), mas a tocasse, nv continuaria sendo 0, porém, como tal condição seria
51
uma situação limite (o mais próxima possı́vel de enlace do ponto), o sistema se-
ria marginalmente estável, apresentando resposta harmônica a entradas degrau. O
enlace corresponderia a uma resposta oscilatória divergente, já que nv = 2 im-
plicaria em nf = 2. A figura abaixo apresenta as diferentes respostas temporais
associadas à localização da curva de Nyquist em relação ao ponto crı́tico, válida
para os casos onde np = 0:
Figura 39: Resposta temporal de sistemas em malha fechada que apresentem np =
0.
Casos de sistemas de controle cuja função de transferência em malha aberta
possui apenas pólos estáveis (np = 0) são recorrentes em uma gama grande de
aplicações tecnológicas. Por esta razão, a previsão de suas respostas temporais por
diagramas como o da Figura 39 é uma boa abordagem de análise. Contudo, existe
uma forma de analisar a estabilidade e o comportamento dessa classe de plantas
de um modo mais prático: através de parâmetros conhecidos como margens de
estabilidade.
Considere uma planta,G(jω), cujo Diagrama de Nyquist seja traçado segundo
a figura a seguir:
52
Figura 40: Curva de Nyquist de planta arbitrária de np = 0.
É desejável que a curva de Nyquist de G nunca toque ou envolva o ponto
crı́tico (−1; 0). Ou seja, se ela nunca passar pela referida coordenada, garante-se
que ela também não a evolverá. Basta impor a seguinte condição pra que o critério
de estabilidade BIBO valha:
G(jω) 6= −1 (82)
O ponto crı́tico possui um módulo de valor unitário, 1, e uma fase rasa, de
−180o. Por isso, é suficiente queG não assuma simultaneamente as duas condições
para que o sistema em malha fechada seja estável:
{¬ [(|G| = 1) ∧ ( G = −180o)]} =⇒ (o sistema é BIBO estável) (83)
Dada a condição de estabilidade, basta verificar se ela vale para a curva de
Nyquist da Figura 40.
De inı́cio, verifica-se a condição de módulo, |G| = 1, respondendo à seguinte
questão: na frequência onde o módulo de G é unitário, sua fase atinge −180o?
53
Figura 41: Análise gráfica de G quando seu módulo é unitário.
Na ilustração acima, percebe-se que o ponto da curva onde |G| = 1 não possui
fase igual a−180o e distaMf , em termos de ângulo, do ponto (−1; 0). Como essa
quantidade expressa uma ”folga” em termos de fase para deixar o sistema longe do
ponto crı́tico, Mf é conhecida por margem de segurança em fase ou simplesmente
margem de fase.
Após analisar a condição de módulo, vale avaliar também a questão do valor
crı́tico da fase de G, respondendo à seguinte pergunta: quando (G) = −180o, o
módulo da planta é unitário?
54
Figura 42: Análise gráfica de G quando sua fase é −180o.
O inverso do módulo da planta quando o fasor a ela associado apresenta
uma orientação de −180o no plano complexo, k, é conhecido como margem
de segurança multiplicativa ou margem de ganho. É o número pelo qual G
ainda pode ser multiplicada antes da sua curva de Nyquist envolver a coorde-
nada (−1; 0). Essa quantidade desempenha o mesmo papel que o ganho assumia
quando se analisava a estabilidade do sistema de controle no Diagrama do Lugar
das Raı́zes.
Em plantas que atendem o critério de estabilidade BIBO, de funções de trans-
ferência próprias (com um número de pólos maior ou igual o número de zeros) e
que detenham todos os zeros no semiplano esquerdo do plano s, se a margem de
ganho é maior ou igual a 1 e, simultaneamente, a margem de fase for semidefinida
positiva, os sistemas de controle em malha fechada a elas associadas serão estáveis
ou marginalmente estáveis. Por outro lado, caso a margem de ganho seja menor
do que 1 e, ao mesmo tempo, a margem de fase for negativa, o sistema de controle
em malha fechada correspondente será instável. Uma G que tem todos os zeros e
pólos localizados no semiplano esquerdo do Diagrama de Argand Gauss e apre-
senta um número de pólos maior ou igual do que o número de zeros é chamada
55
de sistema de fase mı́nima. Tanto a ilustração deste tipo de planta quanto a sua
conexão com a análise de estabilidade do sistema de controle correspondente são
mostrados a seguir:
Figura 43: Mapa de pólos e zeros de um sistema de fase mı́nima.
G é um sistema de fase mı́nima
⇓{
[(log(k) > 0) ∧ (Mf > 0)] =⇒
(
G
1 +G
é BIBO estável
)}
∧
{
[(log(k) < 0) ∧ (Mf < 0)] =⇒
(
G
1 +G
é instável
)} (84)
A proposição 84 estabelece que o logaritmo da margem de ganho e a margem
de fase serem simultaneamente positivos é condição suficiente para que o sistema
de controle em malha fechada apresente todos os pólos no semiplano esquerdo,
ou seja, na região estável do plano s. Ela afirma, ainda, que se essas métricas
forem negativas, pode-se inferir que o sistema em malha fechada é instável. Es-
tas últimas acertivas valem para sistemas de controle lineares que possuam sua
função de transferência de malha aberta de fase mı́nima, mas elas não garantem
nada sobre a estabilidade do sistema em malha fechada caso as referidas quanti-
dades apresentem sinais opostos; neste último caso, deve-se voltar ao diagrama
de Nyquist de G e contar o número de voltas em torno do ponto (−1; 0) no sen-
tido horário para averiguar se existem pólos instáveis na função de transferência
G
1 +G
:
56
log(k)×Mf < 0 =⇒ falta informação para localizar os pólos de
G
1 +G
(85)
Para fins de projeto, a visualização dessas margens através do Diagrama de
Nyquist não é muito prática. A captura de seus valores se torna mais simples
através de um par de gráficos que constitui o que se conhece por diagrama de
Bode.
O diagrama de Bode se desmembra em dois gráficos em escala logarı́tmica: o
gráfico da fase de G e o gráfico do logaritmo do módulo, ambos com o logaritmo
da frequência no papel de abscissa. Nele, é possı́vel vizualizar as duas margens
de forma imediata:
Figura 44: Diagrama de Bode de G.
O diagrama de Bode fornece informações do regime transiente do sistema
de controle em malha fechada: o fator de amortecimento, ζ , e a sua frequência
natural, ωn, cujos valores exatos são extraı́dos do diagrama somente para sistemas
57
de segunda ordem. Porém, uma maneira mais prática de avaliar essas quantidades
é derivar seus valores aproximados, mostrados no diagrama que segue:
Figura 45: Aproximações das métricas do regime transitório, válidas para plantas
com fator de amortecimento ζ 6 0,6.
Para plantas que apresentam fator de amortecimento menor ou igual a 0,6 ,
valem as aproximações explı́citas na Figura 45. A frequência natural do sistema
em malha fechada possui valor próximo ao da frequência de cruzamento do ganho
com 0 dB, enquanto a margem de fase está em torno de cem vezes o fator de
amortecimento:
ωn ≈ ω||G|=1 (86)
58
ζ ≈ Mf
100
(87)
ω||G|=1 é conhecida por banda passante de G, uma vez que a maioria dos sis-
temas de engenharia apresentam uma resposta de fácil medição quando excitados
por sinais de frequência até o referido valor; quando são solicitados por sinais
de frequências acima da banda passante, a resposta,normalmente, possui baixa
amplitude e, caso a atenuação do sinal seja vertiginosa, é de difı́cil medição. Ou
seja, o sistema não atenua muito sinais que entram no sistema com frequências
menores ou iguais às da banda passante.
Além disso, a banda passante de G e do sistema em malha fechada correspon-
dente à mesma planta,
G
1 +G
, dentro das hipóteses já mencionadas quanto às ca-
racterı́sticas fı́sicas de G, são praticamente as mesmas. Como exemplo, suponha
uma G da seguinte forma:
G =
4
s (s+ 2)
(88)
Pela estrutura algébrica de G, é evidente que a frequência natural e o fator
de amortecimento de
G
1 +G
valem ωn = 2 rad/s e ζ = 0,5 , respectivamente.
Plotemos os diagramas de bode de magnitude de ambos os sistemas, o de malha
aberta e o de malha fechada, afim de verificar se as bandas passantes realmente
são próximas:
Figura 46: Banda passante da planta e do sistema em malha fechada.
Olhando novamente a figura acima, mas focando no gráfico de
G
1 +G
, em
59
frequências que não ultrapassam o valor de sua banda passante (na prática, o ωn),
percebe-se que o módulo do sistema em malha fechada,
∣∣∣∣ G1 +G
∣∣∣∣, apresenta valor
unitário (correspondente a 0 dB), enquanto frequências maiores do que ela cor-
respondem a valores de
∣∣∣∣ G1 +G
∣∣∣∣ menores do que a unidade:
Figura 47: Exemplo de coordenada fora da banda passante do sistema em malha
fechada.
Um exemplo de frequência acima da banda passante é do ponto de abcissa
ω = 10 rad/s e ordenada
∣∣∣∣ G(j10)1 +G(j10)
∣∣∣∣ = 10−27,820 ≈ 0,04. Ou seja, se entrarmos
com uma harmônica de frequência 10 rad/s no sistema em malha fechada, sua
saı́da será uma harmônica com amplitude modulada pelo ganho correspondente,
mas multiplicada por 0,04:
G(s)
1 +G(s)
=
y(s)
yc(s)
(89)
60
Figura 48: Relação de amplitudes da entrada e da saı́da na frequência 10 rad/s.
Conforme o esperado, a entrada, yc, de amplitude 100 [unidades], por oscilar
na frequência 10 rad/s, se reduz a praticamente 4% da sua amplitude original
quando é transformada em y pela função de transferência em malha fechada.
3.2.2 Redes proporcional, em avanço de fase e em atraso de fase
Na prática, projetar um sistema de controle no domı́nio da frequência se traduz
na sı́ntese de um compensador, C(s), capaz de moldar o diagrama de bode do sis-
tema em malha aberta, CG(s), para que este apresente banda passante e margem
de fase correspondentes aos requisitos de comportamento em regime transiente
de
CG(s)
1 + CG(s)
. Ao final do projeto, também, o diagrama de bode do módulo
de CG(s) deve apresentar valores de |CG(jω)| em baixa frequência condizentes
com a exigência de erro de rastreamento em regime permanente. Para uma planta
G, o referido sistema de controle é ilustrado na figura a seguir:
Figura 49: Sistema de controle em malha fechada.
61
3.2.2.1 Compensador proporcional
Para o ajuste da banda passante do sistema, basta compensá-lo através de um
ganho K, atribuindo a este último um valor constante, sem a necessidade de
inserção de pólos e/ou zeros na malha aberta. Com isso, apenas a curva do di-
agrama de Bode do módulo se altera, mantendo, assim, a curva da fase inalterada.
Contudo, as margens de estabilidade se alteram com essa mudança no ganho; com
o aumento do ganho, a banda passante se eleva e as margens de ganho e de fase
diminuem, efeitos que podem ser visualizados graficamente no Diagrama de Bode
da figura abaixo:
Figura 50: Efeito do ajuste do ganho de malha aberta na banda passante.
3.2.2.2 Rede em avanço
Quanto à margem de fase requerida, é possı́vel atingı́-la através de um com-
pensador que adicione à fase da função de transferência de malha aberta a fase
que falta para completar o valor. E um compensador indicado para promover
essa mudança é a rede em avanço de fase, cuja estrutura algébrica é mostrada na
sequência:
62
Cavanço(s) =
√
1 + (αTωϕ)
2
1 + (Tωϕ)
2
1 + Ts
1 + αTs
(90)
ωϕ é a frequência de projeto na qual ocorrerá a maior adição de fase, sendo
o ı́ndice ϕ a capacidade máxima do compensador em ângulo, que será definida
pelo projetista mediante um requisito de amortecimento. Além disso, Cavanço é
propositalmente normalizado na frequência ωϕ para que o gráfico da curva de
Bode do módulo da malha aberta permaneça inalterado no ponto do avanço de
fase ϕ. As constantes α e T dependem das referidas variáveis de projeto ϕ e ωϕ
segundo as seguintes relações:
α =
1− sen(ϕ)
1 + sen(ϕ)
(91)
T =
1
ωϕ
√
α
(92)
α possui um valor positivo menor do que 1, o que torna o zero do avanço de
fase mais importante do que o seu pólo (quanto mais próximo o elemento estiver
do eixo imaginário, mais importância ele terá na resposta). E, apesar da rede em
avanço de fase prover amortecimento ao sistema de controle (pelo consequente
aumento da margem de fase), ele pode, dependendo da planta a ser compensada,
aumentar a sensibilidade da malha fechada ao ruı́do. Isso devido à elevação do
ganho do sistema em malha aberta em frequências maiores do que ωϕ, implicando,
assim, no aumento da banda passante do sistema. Tais caracterı́sticas podem ser
observadas na resposta em frequência de Cavanço(jω) e no mapa dos pólos e zeros
da rede, ilustrados na figura abaixo:
Figura 51: Resposta em frequência e mapa de pólos e zeros de uma rede em
avanço de fase.
63
3.2.2.3 Rede em atraso
Finalmente, para atender requisitos de estado estacionário, para que o sistema
de controle atenda determinados nı́veis de erro em regime permanente, o compen-
sador mais adequado a essa tarefa é a rede em atraso de fase, um compensador de
primeira ordem mostrado a seguir:
Catraso(s) =
s+ a1
s+ a2
(93)
a1 e a2 são, respectivamente, as magnitudes do zero e do pólo da rede. Adi-
cionalmente, a1 é maior do que a2, fato esse que torna o pólo mais importante do
que o zero no compensador em atraso de fase. Tal inversão da importância dos e-
lementos da rede implica em dois efeitos derivados da rede em atraso que diferem
do verificados na rede em avanço: no diagrama de Bode da magnitude, a rede
em atraso de fase provê ganhos altos em baixa frequência e ângulos negativos na
curva da fase, conforme pode ser observado na figura a seguir:
Figura 52: Resposta em frequência e mapa de pólos e zeros de uma rede em atraso
de fase.
A subtração ou atraso de fase ocorre de forma apreciável em torno de uma
frequência que se encontra entre as magnitudes do zero e do pólo, especificamente
na média geométrica destes últimos, como indicado na Figura 52. A quantia de
ângulo que é subtraı́da, denotada aqui por $, é calculada com base nos próprios
tamanhos do zero e do pólo segundo a seguinte expressão:
$ = arcsen
(
1− a2
a1
1 + a2
a1
)
(94)
64
O aumento do ganho de malha aberta em baixa frequência, caracterı́stica mais
notável deste tipo de compensação, é o fator determinante para a diminuição do
erro de regime estacionário. Vale a pena, neste momento, equacionar esse erro de
regime permanente para alguns comandos tı́picos em testes de controladores para,
então, esclarecer a sua conexão com a sı́ntese da rede em atraso de fase.
Figura 53: Erro de regime estacionário a comandos tı́picos de teste.
A avaliação do comportamento de grandezas de um sistema LTI em regime
permanente se dá através do cálculo do limite delas quando o tempo t tende ao
65
infinito. Uma maneira equivalente de abordar o mesmo problema, porém, no
plano s, é avaliá-las quando a variável s tende a zero. É dentro deste tipo de
análise que é enunciada a equação do erro de regime permanente de um sistema
de controle de compensador C e de planta G, ambas as funções de transferência
em série no circuito de controle descrito na Figura 49:
e∞ = lim
t→∞
(yc − y) =

1
1+kp
, se yc = 1
1
kv
, se yc = t;
1
ka
, se yc = t
2
2
(95)
As funções que yc assume, yc = 1, yc = t e yc = t2/2, são denotadas, res-
pectivamente, por ”função degrau”, ”sinal rampa”e ”entrada em aceleração” e as
constantes kp, kv e ka são calculadas pelas seguintes equações:
kp = lim
s→0
CG(s) (96)
kv = lim
s→0
sCG(s) (97)
ka = lim
s→0
s2CG(s) (98)
Visto isso, um exemplo simples de projeto de um compensador em atraso de
fase será resolvido, principalmente, para elucidar suas propriedades quanto ao
ganho e à fase e justificar os apontamentos da vantagem e desvantagem que elas
carregam, estando estes últimos destacados na Figura 52.
Seja uma planta G cuja função de transferência é dada a seguir:
G =
100(s+ 1)
s(s+ 20)
(99)
Considere, agora, o projeto de um compensador C que deva atender o seguinte
requisito de regime permanente:
(ka)exigido = 50 (100)
Se o requisito dado é ka, significa que é exigido que o sistema de controle
siga um comando ”em aceleração”, conforme a definição do erro em regime esta-
cionário presente na Equação 95.
66
Inicialmente, vale verificar se a própria planta, G, no papel de função de trans-
ferência de malha aberta é suficiente para seguir a ”entrada em aceleração” com
um erro finito:
C(s) = 1 =⇒ (ka)atual = lims→0 s
2CG(s) = 0× 100
20
= 0 (101)
(e∞)atual =
1
(ka)atual
=
1
0
→∞ (102)
Pelo resultado do cálculo efetuado na Equação 102, percebe-se que um valor
nulo para ka implica em um erro em regime estacionário exorbitante. E para que
este erro assuma um valor finito, é necessário que a função de transferência de
malha aberta, CG(s), possua ao menos mais um s no denominador; só assim,
o s da conta efetuada na Equação 101 será cancelado, permitindo que ka atinja
um valor não-nulo e, consequentemente, que o erro estacionário alcance um valor
bem definido. Portanto, C(s) será redefinida, porém, com um integrador, apenas.
Recalculando, então, o valor de ka segundo o ajuste mencionado, tem-se que:
[C(s)]integrador =
1
s
=⇒ (ka)ajuste = lims→0 s
2G(s) [C(s)]integrador =
100
20
= 5 (103)
Ainda que ka não tenha atingido, até então, o valor do requisito do projeto
(observação oriunda da comparação entre os resultados das Equações 103 e 100),
essa grandeza assumiu um valor não-nulo, garantindo um erro finito. É agora,
então, que entra em cena a rede em atraso de fase. E para a sı́ntese desta última,
faz-se necessário o diagrama de Bode do sistema em malha aberta, cuja função de
transferência, nesta etapa de projeto, se encontra no seginte formato:
G(s) [C(s)]integrador =
100(s+ 1)
s2(s+ 20)
(104)
Esboçando o diagrama de Bode da função de transferência de malha aberta da
Equação 104, chegamos à seguinte resposta em frequência do sistema:
67
Figura 54: Respostas em frequência de G(s) [C(s)]integrador, [C(s)]atraso e{
G(s) [C(s)]integrador
}
× {[C(s)]atraso}.
A resposta em frequência de G(s) [C(s)]integrador da Figura 54, correspondente
às curvas vermelhas do módulo e da fase, mostra que há cruzamento do ganho com
0 dB, mas exibe uma fase que não intercepta a abscissa −180◦. Isso mostra que
não há, portanto, uma margem de ganho finita para se medir no referido sistema,
essa margem é ”infinita”: não importa o quanto você aumente o ganho de malha
aberta, G nunca atingirá o ponto (−1; 0). A margem de fase, entretanto, é finita
e possui um valor desejável (maior do que 60◦ como boa prática de projeto), de
aproximadamente 65◦.
O projeto do compensador em atraso de fase que atuará sobreG(s) [C(s)]integrador
visa atender o ka requerido e, ao mesmo tempo, garantir que sua compensação não
deteriore o amortecimento do sistema, que não diminua, portanto, a sua margem
de fase, preocupação essa decorrente do fato deste compensador causar subtração
68
de ângulo na malha aberta.
Um decaimento grande da margem de fase ocorrerá no processo somente se
a frequência de cruzamento com 0 dB estiver na faixa espectral de operação
desse compensador. E para contornar essa situação, basta escolher essa faixa de
operação, [a2; a1] (parâmetros presentes na Equação 93), para que ela não con-
tenha a frequência na qual se mede a margem de fase, mas, ao mesmo tempo,
em regiões de baixa frequência para atender ka. A curva amarela, da Figura 54,
mostra uma boa escolha para a faixa de operação da rede em atraso afim de aten-
der as demandas aqui mencionadas. Ou seja, basta impor um valor a a1 muito
menor do que a frequência sobre a qual a margem de fase está definida.
Como a frequência de cruzamento com 0 dB está em torno de 5 rad/s, o
tamanho do zero do compensador, a1, para que seja muito menor que a referida
frequência, será tomado como o valor dez vezes menor que o dela:
a1 =
5
10
(105)
E para criar um vı́nculo entre o tamanho do pólo do compensador, denotado
por a2, e o valor do requisito de regime estacionário, (ka)exigido, inserimos a rede
em atraso de fase na malha e reescrevemos a equação de ka com o novo sistema
em malha aberta:
(ka)exigido = lims→0
s2
{
G(s) [C(s)]integrador
}
× {[C(s)]atraso} (106)
Inserindo os valores de a1 e (ka)exigido, presentes nas Equações 105 e 100, res-
pectivamente, na Equação 93, que contém a estrutura algébrica de [C(s)]atraso, e
alocando o resultado na Equação 106, estando esta última com a funçãoG(s) [C(s)]integrador
definida da Equação 104, o valor de a2 pode ser calculado, chegando-se ao seguinte
valor:
a2 =
5
100
(107)
A rede em atraso, então, toma a seguinte forma:
[C(s)]atraso =
s+ 5
10
s+ 5
100
(108)
A resposta em frequência do sistema compensado após a inserção da rede em
atraso, representada pelas curvas em rosa de
{
G(s) [C(s)]integrador
}
×{[C(s)]atraso}
na Figura 54, mostra que houve um pequeno abaixamento da margem de fase
69
(que foi para aproximadamente 60◦), porém, o decaimento da fase em baixas
frequências ocasionou o cruzamento da curva da fase com a abscissa −180◦.
Nessa interceptação, a margem de ganho, k, é visivelmente menor do que 1. E esta
margem de fase positiva aliada ao logaritmo da margem de ganho negativo (neste
caso, −23 dB) fazem recair sobre este sistema a condição 85: não se pode dizer
nada sobre a estabilidade do sistema em malha fechada, não é possı́vel inferir se
os pólos de
{
G(s) [C(s)]integrador
}
× {[C(s)]atraso}
1 +
{
G(s) [C(s)]integrador
}
× {[C(s)]atraso}
estão ou não localizados
na região instável do plano s.
Então, para verificar a estabilidade desse sistema de controle, é necessário
traçar o Diagrama de Nyquist de
{
G(jω) [C(jω)]integrador
}
× {[C(jω)]atraso}, exi-
bido na figura que segue:
Figura 55: Diagrama de Nyquist de
{
G(jω) [C(jω)]integrador
}
× {[C(jω)]atraso}.
Na Figura 55, o ponto (−1; 0) é enlaçado duas vezes, uma vez no sentido
horário e outra no sentido anti-horário. E as voltas que a referida curva faz ao
70
redor da coordenada em questão são destacadas a seguir:
Figura 56: Envolvimentos da coordenada (−1; 0) pela curva de Nyquist.
Da Figura 56, contabilizando as duas voltas, temos a resultante dos enlaces
em torno do ponto crı́tico: nv = (+1) + (−1) = 0. Além disso, o número de
pólos de
{
G(jω) [C(jω)]integrador
}
×{[C(jω)]atraso} localizados na região instável
do plano s é nulo: np = 0. Portanto, aplicando o Critério de Nyquist ao problema,
constata-se que a soma de nv com np é nula, valor correspondente ao número de
pólos instáveis do sistema em malha fechada:
nf = np + nv = 0 =⇒
{
G(s) [C(s)]integrador
}
× {[C(s)]atraso}
1 +
{
G(s) [C(s)]integrador
}
× {[C(s)]atraso}
é estável
(109)
A resposta temporal desse sistema em malha fechada é mostrada a seguir:
71
Figura 57: Sistema de controle rastreando um comando ”em aceleração”.
Depois de averiguar que o sistema em malha fechada é estável, que ele atende
o requisito de regime permanente e que, além disso, também exibe uma boa res-
posta transitória, chega-se ao fim do projeto do compensador C para a planta G
desse exemplo:
C(s) =
{
[C(s)]integrador
}
× {[C(s)]atraso} =
1
s
(
s+ 5
10
s+ 5
100
)
(110)
3.2.3 Altitude e velocidade
As dinâmicas restantes