Baixe o app para aproveitar ainda mais
Prévia do material em texto
UNIVERSIDADE FEDERAL DE UBERLAˆNDIA APOSTILA MECAˆNICA DO VOˆO Controle Automa´tico de Aeronaves Autor: Alexandre MASSON Supervisor: Dr. Leonardo SANCHES 2017 2 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 comple- mentares a esta apostila que podem ser combinados com esta u´ltima 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 quantidades cinema´ticas sera˜o medidas e decompostas (em 3 direc¸o˜es), para avaliar o movimento da aeron- ave 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ω2L −CueD CueL −(aL + bLα) −(aD + bDα) −Cω2L Cω2D −CueL CueD ( cosα senα ) (29) A` excessa˜o de Va, α e β, que variam no tempo, os demais termos sa˜o con- stantes 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 req- uisitado 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 asso- ciados 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 semovimentar 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 1V ∗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 si 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, 5pi 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 da diag- onal. 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 denom- inador, 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[pi−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 nen- hum valor particular, tampouco diverge: 23 Figura 11: Eixo imagina´rio como limite de estabilidade dosistema. 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) = pi180rad, foi plotada com a ajuda do MATLAB: Figura 13: Resposta de δφ(s) δφc(s) a uma entrada degrau de pi 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) = pi 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 cont- role. 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´losde 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 valores de 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 de- sprezado, 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,2pi 180 sen(2pi100t) 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−i2piftnβ(t)dt = 0,2 [ i2pif (i2pif)2 + (2pi100)2 ] (71) δβc(f) = ∫ ∞ 0 e−i2piftδβc(t)dt = e−4(i2pif) − 1 4(i2pif)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ˆnica constitui 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 lin- eares 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 fechadae´ 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˜o 77 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 re- stringe 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 partic- ular: 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 ap- resenta 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; nesteu´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 car- acterı´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 corre- spondem 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 el- ementos 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, re- spectivamente, 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, re- spectivamente, 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
Compartilhar