Identificação de Sistemas Dinâmicos Lineares

IDENTIFICAÇÃO DE SISTEMAS DINÂMICOS LINEARES ANTÔNIO AUGUSTO RODRIGUES COELHO LEANDRO DOS SANTOS COELHO processo dispositivo de medição modelo parâmetros IDENTIFICAÇÃO DE SISTEMAS DINÂMICOS LINEARES UNIVERSIDADE FEDERAL DE SANTA CATARINA Reitor Lúcio José Botelho ViceReitor Ariovaldo Bolzan EDITORA DA UFSC Diretor Executivo Alcides Buss Conselho Editorial Eunice Sueli Nodari Presidente José Isaac Pilati Luiz Henrique de Araújo Dutra Luiz Teixeira do Vale Pereira Sérgio Fernando Torres de Freitas Tânia Regina Oliveira Ramos Vera Lúcia Bazzo Antonio Augusto Rodrigues Coelho Leandro dos Santos Coelho IDENTIFICAÇÃO DE SISTEMAS DINÂMICOS LINEARES Editora da UFSC Florianópolis 2004 dos autores SUMÁRIO SUMÁRIO 44 Modelagem de Processos via Resposta em Frequência 83 PREFÁCIO Controle Contínuos e Discretos As simulações através de computador são indispensáveis para facilitar o entendimento dos diferentes algoritmos de identificação de processos Muitas das simulações apresentadas foram realizadas com o auxílio do software científico denominado Matlab mas outros softwares aplicativos em controle com características similares também podem ser utilizados Sistema Modelagem e identificação Descrição de sistemas modelos contínuos Descrição de sistemas modelos discretos Problemas 11 SISTEMA O conceito de sistema pode ser definido de diferentes formas Em controle de processos denotase como um objeto ou uma coleção de objetos que realiza um certo objetivo e cujas propriedades pretendese estudar Alguns exemplos são sistema de fabricação de papel ou cerâmica sistema solar circuito elétrico servomecanismo de posição sistema biológico ou econômico manipulador robótico reator coluna de destilação entre outros Dorf Bishop 1995 Ogata 2002 A figura 11 ilustra os principais elementos de um sistema de controle entrada causa u h saída efeito y H Figura 11 Componentes de um sistema de controle Os problemas de controle associados à estrutura da figura 11 são Análise é conhecida a entrada u o sistema h e devese obter a saída y Estes procedimentos propiciam a obtenção de modelos que representam a dinâmica do sistema processo ou planta Para fins de controle de processos não se pretende encontrar um modelo matemático exato mas um modelo adequado para uma determinada aplicação Na prática utilizase a hipótese básica para a elaboração de modelos de que processos reais em geral não necessitam obrigatoriamente de modelos complexos Hang e Chin 1991 Gessing 1996 Ljung 1999 O modelo de um sistema é uma equação matemática utilizada para responder a questões sobre o sistema sem a realização de experimentações através de um modelo podese calcular ou decidir como o sistema comportase sob determinadas condições operacionais A utilização do modelo para simulação do sistema constituise um procedimento de baixo custo e seguro para experimentar o sistema Entretanto a validade adequação dos resultados de simulação depende completamente da qualidade do modelo matemático do sistema Alguns dos diferentes propósitos para a utilização de modelos matemáticos em automação industrial são Previsão é uma tentativa de prever os estados futuros do sistema comportamento dinâmico e está limitada à precisão do modelo e aos efeitos das perturbações atuantes presentes no sistema Análise e projeto de sistemas de controle proporciona um vasto campo para aplicação em modelagem e identificação na sintonia de controladores clássicos síntese de algoritmos de controle adaptativos e preditivos e na estimação do estado de variáveis nãomensuráveis por exemplo a estimação da velocidade a partir da posição e uma medida indireta Supervisão utiliza a simulação com base no modelo matemático para a avaliação das características operacionais do sistema para o projeto de engenharia ou para o treinamento de operadores Muitas vezes é também utilizado na detecção de erros e diagnóstico Otimização empregado na tomada de decisões nos mais variados campos no escalonamento na manutenção e na economia em sistemas industriais maximizar produção minimizar custos etc A otimização de sistemas necessita de modelos matemáticos precisos A resposta impulsiva está relacionada com a função de transferência por ht L1Hs 14 onde L1 é a transformada inversa de Laplace e t é o tempo contínuo A função de transferência está relacionada com a representação de estados no caso monovariável de acordo com a seguinte relação Hs csI A1b 15 Com base na lei das tensões de Kirchhoff temse a equação que descreve a dinâmica do sistema elétrico ou seja Capítulo 1 Introdução Tabela 11 Código em Matlab do sistema mecânico da equação 19 A relação WsVs é dada por Ws KaJsBLsRKa Kb A representação matemática e o comportamento dinâmico do motor DC são conhecidos se os parâmetros são fornecidos pelo fabricante especificações de catálogo iii Representação matemática do sistema monotanque Considere o sistema monotanque conforme a figura 17 onde A é a área do tanque m² a é a área do tubo de saída m² ht é o nível do líquido no tanque m ut é a vazão de entrada m³seg e qt é a vazão de saída m³seg Para avaliar suas características operacionais devese obter uma equação diferencial nãolinear para o sistema calcular um modelo matemático relacionando as variáveis de entrada ut e saída ht Assim utilizase a lei de Bernoulli que descreve a relação entre a velocidade de vazão da saída mseg e o nível do líquido no tanque isto é vt 2ght onde g é a aceleração da gravidade A equação relacionando o vazão da saída qt e a velocidade de vazão da saída vt é dada por qt avt O volume de líquido no tanque em um instante t é calculado por Aht m³ e modificase de acordo com a diferença entre o fluxo de entrada e saída denominado balanço de massa isto é dAhtdt qt ut Substituindose as equações 112 e 113 na equação 115 obtémse uma equação diferencial nãolinear para avaliar o comportamento do nível do líquido ou seja dhtdt a2gAht 1A ut Pelo conhecimento da dimensão dos diferentes elementos que compõem o sistema monotanque podese avaliar sua dinâmica pela equação 116 por exemplo é possível determinar o nível ht quando o fluxo de entrada ut é conhecido O fluxo de saída é calculado por qt a2ght A seguir apresentamse de forma resumida as três abordagens utilizadas para descrever sistemas discretos Com o auxílio do software científico Matlab podese elaborar e converter os modelos contínuos em discretos Também é possível através das equações que descrevem as transformações retangular ou trapezoidal determinar os correspondentes modelos contínuos e discretos i Função de transferência discreta A função de transferência discreta é a relação entre a transformadaz da saída Yz pela transformadaz da entrada Uz isto é Hz YzUz A relação Hz é uma razão de dois polinômios em z e está representada por Hz BdzAdz onde Bdz bj zj para j0 até m Adz ai zi para i0 até n sendo n m adn l e n é a ordem do sistema Os elementos bj adi n com j 0m e i 0n são desconhecidos e devem ser determinados através de uma modelagem matemática do sistema ou através da identificação Algumas vezes assumese que o valor de n é conhecido a priori ii Resposta impulsiva discreta A resposta impulsiva está relacionada com a função de transferência por ht Z1Hz onde Z1 é a transformadaz inversa e t é o tempo discreto 142 EXEMPLO DE MODELAGEM POR ANÁLISE EXPERIMENTAL i Processo discreto de primeira ordem Considere o processo caracterizado pela seguinte equação a diferenças yt 1 θ yt ut 128 onde θ é o parâmetro desconhecido e yt 0 t 0 Admita o seguinte modelo para estimação ỹt 1 ŷθt ut 129 sendo ŷ a estimativa de θ e ỹt 1 a previsão da saída ou o valor da saída no instante t 1 baseado no parâmetro estimado ŷ Seja a função custo dos mínimos quadrados dada por Jt 12 Σ e² k 130 onde et yt ỹt 131 Substituindose as equações 128 e 129 na equação 131 obtêmse o erro de estimação em função das medidas de entrada e saída atual e anteriores do processo de acordo com et θ yt 1 ỹt 1 yt ut 1 ỹ yt 1 132 Supor para um certo sistema que não é conhecido o parâmetro θ mas as medidas de entrada e saída no intervalo de tempo da experimentação 0 t N Assim obtémse o estimador dos mínimos quadrados pela diferenciação de Jt em relação a θ resultando O tempo discreto t é utilizado para indicar que a estimação está baseada nas medidas até e incluindo o instante t A estimação do parâmetro θ dado pela equação 133 minimiza 130 sob a suposição que o processo está representado pela equação 128 De acordo com a equação 133 observase que é um algoritmo nãorecursivo estima somente um parâmetro e necessita de todas as medidas de entrada e saída da experimentação procedimento offline A solução da equação 133 estimativa de θ pode ser obtida com um software numérico tal como o Matlab A tabela 12 ilustra a programação em Matlab da equação 133 onde se admite θ 09 uma sequência grau unitário e o tempo de experimentação com 30 amostras 15 PROBLEMAS 1 Considerando a constante de tempo mecânica τm JB maior que a constante de tempo elétrica τe LR do motor DC controlado por armadura obter a equação diferencial de primeira ordem em termos da velocidade e tensão de armadura 2 De acordo com as especificações de um fabricante que comercializa um motor DC temse os seguintes parâmetros R 20 Ω J 2 Nmseg²rad B 01 Nmsegrad Ka 1 NmA Kb 3 Vsegrad Avaliar o comportamento do processo por simulação para as seguintes condições i Admitir condições iniciais nulas τm τe e aplicar um grau unitário na tensão de armadura para observar o comportamento dinâmico da velocidade ii Observar as dinâmicas temporais analógica e digital período de amostragem de 01 seg sob as mesmas condições de operação do item i 3 Avaliar por simulação o comportamento do sistema de nível equação 116 quando ut 1 t 0 grau unitário e h0 0 Observar também a característica temporal para ut 0 e h0 025 m Admitir A 1 e a 2g 1 Adicionalmente considerar o ponto de operação u0 h0 e mostrar que o modelo linearizado é de primeira ordem relacionando HsUs 4 Simular os comportamentos de carregamento e descarregamento da tensão no capacitor do circuito RC a partir dos modelos contínuo e discreto Admitir R 1 MΩ C 1 μF e Ts 01 e 1 seg Seja a tensão de entrada do circuito de acordo com ut 1 0 t 10 0 fora 5 Implementar o programa em Matlab da tabela 12 para estimativa do parâmetro θ relativo ao processo discreto de primeira ordem representado por yt θyt1 ut1 onde θ 075 ut é uma sequência grau unitário discreto e sendo 20 amostras o tempo total da experimentação 6 A tabela mostra os valores medidos da resposta impulsiva de um sistema Calcular e representar graficamente a resposta ao grau unitário t 0 01 02 03 04 05 06 07 08 09 1 ht 0 95 59 014 27 27 14 03 012 012 001 Para fins de aplicação de controle digital direto qual deve ser o período de amostragem como critérios de informação Bayesiana de Akaike ou minimum description length combinam a variância residual e a ordem do modelo Haber e Unbehauen 1990 A meta do algoritmo de otimização é a minimização de um critério de desempenho Se todas as restrições e condições forem atendidas o modelo encontrado pode ser aceito Caso contrário se uma das condições impostas é violada todos os procedimentos de identificação estimativa de parâmetros e diagnóstico do modelo devem ser reavaliados até que um modelo apropriado seja encontrado Ljung 1999 A área de identificação tem tido considerável interesse nos últimos anos para fins de previsão supervisão diagnóstico e controle Observase a aplicação em diversos campos da engenharia tais como processos químicos sistemas biomédicos sistemas socioeconômicos sistemas elétricos e outros mais Em uma tarefa de identificação diferentes procedimentos para geração do sinal de entrada medição da saída e armazenamento dos dados são utilizados entre os quais Identificação de um processo pelo teste de resposta ao degrau o processo é submetido a uma mudança na entrada do tipo degrau e o armazenamento das medidas é realizado por algum registrador Com a curva de reação do processo é possível aplicar diferentes técnicas gráficas numéricas ou computacionais para modelar o sistema controlado por funções de transferência de primeira ou segunda ordem processos com modelos de ordem reduzida conforme ilustra a figura 22 O teste de resposta ao degrau só tem validade para processos lineares ou para processos não lineares que sejam aproximadamente lineares na vizinhança do ponto de operação O teste de resposta ao degrau por natureza não permite a estimação de modelos de ordem superior já que o sinal degrau tem uma pobre composição em frequência Identificação pelo teste da resposta em frequência o processo é submetido a uma entrada harmônica sinal senoidal De acordo com as curvas de magnitude e fase é possível identificar as frequências de corte avaliandose a influência dos zeros e pólos e consequentemente a função de transferência estimada correspondente conforme mostra a figura 23 Identificação offline com o auxílio de sinais de teste apropriados de entrada ruído branco sequência binária pseudoaleatória excitase o processo e armazenamse as medidas de entrada e saída para aplicação e avaliação a posteriori dos algoritmos de estimação nãorecursivos conforme ilustra a figura 24 Para isso é necessário que a estrutura do modelo seja disponível ou seja é preciso selecionar a ordem do modelo o número de pólos e zeros e eventualmente o valor do atraso de transporte ou tempo morto Todo cálculo dos parâmetros é feito em um momento e o tempo de avaliação é diferente daquele em que se realiza o ensaio Por isso se diz que a estimação é feita offline Para este tipo de abordagem utilizamse na prática modelos discretos para os processos Isso se justifica pelo fato de que os algoritmos de identificação trabalham com os valores de amostras dos sinais de entrada e saída A identificação de modo geral consiste em três etapas determinação da estrutura estimação dos parâmetros e validação do modelo A figura 26 ilustra as diferentes etapas da identificação de processos ZN Y U Y ytt 1 N U utt 1 N onde Y é o conjunto dos sinais de saída medidos yt é o conjunto das entradas do sistema ut N é o universo da experimentação e t especifica os instantes de amostragem múltiplos de Ts Ts 2Ts 3Ts O procedimento da figura 27 pode ser utilizado para identificação offline ou online se os sinais de entrada e saída são amostrados e o cálculo para ajuste do modelo é realizado não recursivamente ou recursivamente respectivamente Uma grande variedade de métodos está disponível na literatura para identificar processos ambos online e offline As técnicas de estimação de parâmetros de processos podem estar apoiadas em sistemas de conhecimento inteligência computacional em tarefas de supervisão auxiliar o usuário na utilização e regulação dos algoritmos e em modelagem Alguns softwares disponíveis comercialmente para modelagem de processos monovariáveis e multivariáveis são ADAPTx AIDA CONTROL STATION CONTSID ICAI MATLABSIMULINK PITOPS SIMNON WINFACT A tomada de decisões e a resolução de problemas são dependentes do acesso à informação adequada sobre o problema a ser resolvido Frequentemente a informação avaliada está originalmente na forma de dados ou observações que requerem interpretação antes da análise e tomada de decisões A derivação de uma relevante descrição do sistema a partir de dados observados é denominada identificação de sistemas processos e a descrição do sistema resultante é denominada modelo Todavia existem diversas formulações para o que é identificação de processos A seguir mencionamse algumas definições estabelecidas na literatura Identificação é a determinação baseada nos sinais de entrada e saída de um sistema em que para uma classe especificada de sistemas o sistema sob teste é equivalente Zadeh 1965 Identificação de sistemas representa a interface entre modelos do mundo matemático e o mundo real Johansson 1993 i 1960 desenvolvimento das raízes estatísticas esse período essencialmente começa com K F Gauss em 1809 e termina por volta de 1960 quando modelos paramétricos explícitos começam a ocasionar interesse na comunidade de controle Durante esse período todos os conceitos estatísticos essenciais usados em sistemas de identificação são desenvolvidos Entre os desenvolvimentos desse período podese citar regressão linear mínimos quadrados método da máxima verossimilhança aproximação estocástica análise por série temporal e método da variável instrumental ii 19601970 proliferação dos métodos de identificação aparecimento acentuado de diversos métodos devido principalmente aos seguintes fatores A equação a diferenças linear básica para as relações entradasaída pode ser escrita como regressão linear e assim o método dos mínimos quadrados pode ser aplicado Desenvolvimento de técnicas espectrais e de correlação aplicadas a séries temporais Desenvolvimento de diversos paradigmas de convolução e expansão de funções de resposta ao impulso iii 19701985 consolidação da teoria e prática em sistemas de identificação as linhas mestras podem ser distinguidas em conexões entre diferentes metodologias e intenso desenvolvimento de software para o tratamento de problemas em identificação iv 1985 surgimento de novas ideias em raízes estatísticas em meados dos anos 80 a visão estatística de sistemas de identificação já está madura e sedimentada As estruturas clássicas de estimação de parâmetros alcançam relativo sucesso e podem ser coerentemente transferidas para o mundo dos sistemas dinâmicos O desenvolvimento de ferramentas computacionais poderosas e o aparecimento maciço de pacotes de softwares comerciais também contribuiram para o aprimoramento da área de identificação Entre os tópicos com desenvolvimentos significativos hoje em dia que na maioria das vezes possuem pequena interação com métodos estatísticos citase métodos de subespaco para modelos espaçoestados identificação para controle propriedades de rejeição e tratamento de ruído modelos blackbox wavelet etc 25 APLICAÇÃO EM CONTROLE ADAPTATIVO Um exemplo de utilização da identificação de sistemas em engenharia de controle de processos que tem sido enfatizado nos últimos anos é na implementação de controladores adaptativospreditivos Justificase a utilização destes algoritmos de controle pela adequação para tratar processos variantes no tempo com parâmetros desconhecidos e adicionalmente em aplicações práticas devido à ineficiência do desempenho dinâmico de algoritmos de controle com ganhos fixos na presença de complexidades dinâmica assimétrica parâmetros variantes nãolinearidades no processo controlado Roffel et al 1989 Malik et al 1991 Gessing 1996 Os controladores adaptativos do tipo autoajustável selftuning apresentam uma malha adicional que automaticamente modifica os parâmetros para compensar as variações ocorridas no processo ou no meio ambiente trazendo benefícios significativos em relação a outros tipos de controladores pois as variações da dinâmica são acompanhadas pelo identificador e consequentemente pelo controlador a cada período de amostragem Entre os controladores adaptativos o regulador autoajustável de variância mínima proposto por K J Åström e B Wittenmark em 1973 tem tido crescente utilização implementação e variação de projeto no contexto de controle de processos Coelho et al 1988 Åström e Wittenmark 1995 A partir desta estrutura de controle foi possível a elaboração de controladores adaptativos robustos isto é o controlador autoajustável de variância mínima generalizada Clarke e Gawthrop 1975 e posteriormente o controlador preditivo generalizado Clarke et al 1987 Esta estratégia de controle está dividida em duas etapas identificação e controle conforme ilustra a figura 210 A aplicação em tempo real dos algoritmos de identificação é interessante para vários propósitos entre os quais rastreamento de parâmetros variantes no tempo detecção diagnóstico filtragem controle adaptativopreditivo e redes neurais Capítulo 2 Noções básicas sobre identificação 39 entrada sistema saída identificador modelo do sistema projeto do controlador parâmetros do controlador implementação do controlador controlador adaptativo Figura 210 Controlador adaptativo do tipo autoajustável A etapa de identificação consiste na obtenção dos parâmetros do modelo do sistema controlado através de um estimador de parâmetros recursivo A segunda etapa de controle consiste no cálculo dos parâmetros do controlador pela minimização de um critério de desempenho especificado procedimento ótimo ou pela técnica de alocação de pólos procedimento clássico Seborg et al 1986 26 PROBLEMAS 1 Ilustrar uma experimentação e descrever a aplicação das técnicas de modelagem baseadas nas respostas ao degrau e frequência identificação offline e online 2 Avaliar a qualidade do modelo estimado em relação ao teste de resposta ao degrau em malha aberta e ao teste frequencial pelo diagrama de Nyquist considerando processo real Gps 1020s 4080 s5 24s4 218s3 956s2 2456s 4080 modelo estimado ˆGps 734 s2 168s 734 3 Avaliar a qualidade dos modelos estimados em relação aos testes das respostas ao degrau em malha aberta e frequencial diagrama de Nyquist para os seguintes processos real e estimados processo real Gps 1 062s 1 7 172s 1 modelos estimados ˆGp1s e243s 198s 1 171s 1 ˆGp2s 1 48s 1167s 1 Adicionalmente validar os modelos estimados com base no comportamento temporal em malha fechada para uma entrada degrau unitário utilizando um controlador PI com duas sintonias isto é ParâmetrosCaso a b c Kc 28 56 663 Ti 17 17 10 4 Uma malha de controle com controlador PI uma válvula nãolinear e um processo de terceira ordem estão ilustrados na figura 211 Figura 211 Malha de controle analógica com controlador planta e válvula A característica da válvula está representada por v fu u4 u 0 Capítulo 2 Noções básicas sobre identificação 41 Verificar por simulação Simulink a necessidade de ressintonização do controlador PI ajustado com Kc 015 e Ti 1 assumindo os seguintes níveis de operação i Mudança de referência de 02 para 03 ii Mudança de referência de 1 para 11 iii Mudança de referência de 5 para 52 A sintonia PI deve garantir uma dinâmica adequada em malha fechada sem sobresinal e com tempo de estabilização de aproximadamente 10 seg Idealizandose uma concepção de controle do tipo PI ganho escalonado gain scheduling quais os ganhos nas três faixas de operação CAPÍTULO 3 MODELOS DE PROCESSOS DE ORDEM REDUZIDA E COMPLEXOS Introdução Modelos contínuos de ordem reduzida Modelos discretos de ordem reduzida Modelos de sistemas discretos complexos Perturbações nos modelos discretos Problemas 31 INTRODUÇÃO Vários processos industriais são bem representados por funções de transferência modelos matemáticos de primeira ou segunda ordem Os sinais podem ser representados utilizandose variáveis contínuas ou discretas e as funções de transferência através das transformadas de Laplace e z respectivamente Seborg et al 1989 Ogata 2002 32 MODELOS CONTÍNUOS DE ORDEM REDUZIDA A representação usual de um modelo matemático para uma planta industrial é a função de transferência de primeira ordem com atraso de transporte isto é Gpls YsUs Keθsτs 1 31 onde s é o operador Laplace s ddt K é o ganho θ é o atraso de transporte contínuo τ é a constante de tempo e Ys e Us são as transformadas de Laplace da saída e da entrada respectivamente A seguir discutese como os parâmetros K θ τ que caracterizam o comportamento dinâmico de um sistema de primeira ordem podem ser identificados com base na curva resposta de reação do processo i Ganho de regime O ganho de um sistema é definido por K ΔytΔut 32 onde Δyt é a variação da saída e Δut é a variação da entrada O exemplo 31 mostra o procedimento de determinação do parâmetro K para um processo de nível Exemplo 31 Considerase o processo de nível de líquido em regime estável conforme a figura 31 onde nas condições operacionais têmse as vazões de entrada e saída constantes e iguais Figura 31 Sistema de nível de líquido As figuras 32 e 33 ilustram os comportamentos do nível e o cálculo do ganho considerando um aumento súbito na forma de degrau na vazão de entrada e um aumento em degrau na vazão de saída através da abertura da válvula de saída respectivamente Observações De acordo com os resultados das figuras 32 e 33 observase que o nível muda quando se realiza uma mudança no processo O ganho de um processo tem unidades específicas que devem ser avaliadas cuidadosamente para cada aplicação particular ii Constante de tempo A constante de tempo para processos industriais caracterizada em uma função matemática de primeira ordem é definida pelo tempo requerido medido a partir do ponto onde o sinal de saída começa a mudar para que a saída do processo de primeira ordem atinja 632 do valor da variação total depois que ocorreu uma mudança na entrada Mollenkamp 1988 Exemplo 32 Dois fluidos quente e frio são misturados num sistema de bombeamento conforme ilustra a figura 34 A temperatura das vazões combinadas é medida em um ponto à jusante do ponto de mistura Qual é a resposta da temperatura medida após um aumento em grau na vazão do fluido frio quente frio vazão de fluido frio temperatura no ponto de mistura temperatura medida θ Figura 34 Mistura de dois fluidos quente e frio Adicionalmente alguns processos industriais são inerentemente bem representados por uma função de transferência de segunda ordem de acordo com as equações 35 ou 36 isto é Gp2s Keθsτ1s 1τ2s 1 35 Gp2s Keθss2ωn2 2ξsωn 1 36 onde τ1 e τ2 são as constantes de tempo K é o ganho θ é o atraso de transporte contínuo ωn é a frequência natural e ξ é o fator de amortecimento A frequência natural determina o período de oscilação da saída e o comportamento dinâmico está caracterizado por ξ ou seja Se ξ 1 sistema sobreamortecido raízes reais e desiguais τ1 τ2 Se ξ 1 sistema criticamente amortecido raízes reais e iguais τ1 τ2 Se 0 ξ 1 sistema subamortecido raízes complexas conjugadas Observações Seja a função de transferência YsUs NsDs As raízes do polinômio do denominador D são denominadas pólos do sistema enquanto que as raízes do polinômio do numerador N são denominadas zeros do sistema Considere o sistema de terceira ordem caracterizado pela seguinte função de transferência Gs w2ns2 2ξwns w2nYs 1 Se a condição 1Ys 0 cw é assegurada então a resposta do sistema de terceira ordem pode ser aproximada pelas raízes dominantes do sistema de segunda ordem Dorf e Bishop 1995 33 MODELOS DISCRETOS DE ORDEM REDUZIDA Embora muitos processos sejam contínuos por natureza os modernos sistemas de controle utilizados para controlar processos baseiamse em computadores digitais e aplicam algoritmos de controle digital Jacquot 1995 O sinal de controle é implementado em intervalos de tempo com o intervalo denominado período de amostragem Ts Se a técnica de projeto de controle por computador baseiase no modelo matemático do processo controlado então um modelo discreto da planta é necessário para projetar e calcular a saída do sinal de controle As diferentes formas de representar uma malha de controle são ilustradas na figura 35 a sistema de controle contínuo analógico b sistema de controle digital por computador c sistema de controle discreto amostrado Figura 35 Concepções de malhas de controle analógico e digital A representação função de transferência discreta equivalente de modelos com funções contínuas de primeira e de segunda ordem pode ser calculada com as técnicas de sistemas amostrados O procedimento de cálculo utiliza o segurador de ordem zero em cascata como a função de transferência contínua do processo e um conjunto de amostras da saída do processo em períodos de amostragem está disponível para identificação e controle conforme ilustra a figura 36 Us Gs h Us Ys Ys Figura 36 Segurador de ordem zero em série com o processo onde Gps é a função de transferência do processo Ghs é a função de transferência do segurador de ordem zero Ys e Ys são os sinais de saída contínuo e pulsado discreto e Us e Us são os sinais de entrada contínuo e pulsado discreto A função de transferência amostrada do sistema é calculada por Yz ZGhsGps1z1 Gpss 37 e sendo Ghs da forma Ghs 1esTss Gpz fracb01 a1 z1 b0 K2 fracK1τ 0133 a1 eaTs 05866 ii Caso com o segurador de ordem zero Gz fracb0 z11 a1 z1 b0 1 eaTs K2τ 04134 a1 eaTs 05866 Podese observar que Grz Grz e assim os casos i e ii conduzem a resultados onde se tem modelos matemáticos diferentes na estrutura e nos parâmetros 34 MODELOS DE SISTEMAS DISCRETOS COMPLEXOS Se um computador é utilizado para controlar um processo o sistema de controle pode ser representado pela configuração apresentada na figura 37 AD conversor analógicodigital DA conversor digitalanalógico Figura 37 Sistema de controle por computador O modelo entradasaída para o processo utilizando a transformada z pode ser obtido como Yz GpzUz Gplz zdb0za1 zd1b01a1z1 onde Gpz é a função de transferência do processo incluindo o segurador de ordem zero e pode ser representada por uma razão de dois polinômios em z1 isto é Gpz zd Bz1 Yz Uz ou Az1yt zdBz1ut A equação 314 é uma representação usual na literatura de identificação e controle digital Os polinômios Az1 e Bz1 relacionados aos pólos e zeros de malha aberta do processo respectivamente são da forma Az1 1 a1 z1 a2 z2 an zn Bz1 b0 b1 z1 b2 z2 bn zn A correspondente equação de diferenças da planta é yt a1 yt1 an ytna b0 utd b1 utd1 bn utdnb onde yt e ut são os sinais de saída e controle respectivamente disponíveis nos instantes de amostragem t nTs onde n ℤ Exemplo 34 Obter a equação a diferenças equivalente e avaliar as características de malha aberta do modelo discreto relativo ao processo de segunda ordem representado por Gps 1 3s 15s 1 dTs θ d1Ts Gp2z zdb0zb1z2a1za2 zd1b0b1z11a1z1a2z2 Converte uma função de transferência continua para discreta Gpsnumsdens Gpznumzdenz Agora admitindose que vt é da forma vt Cz1et Az1 a equação 315 pode ser reescrita como yt Bz1utd Cz1et Az1 316 num1den15 8 1Ts1 Figura 39 Ruído branco com média nula e variância unitária O resultado da figura 39 pode ser obtido executandose o programa em Matlab da tabela 32 numzdenzc2dmnumdenTszoh dispzeros do sistema discreto O modelo matemático CÁRIMA em substituição ao modelo CARMA da equação 316 está caracterizado pela seguinte equação de diferenças Az1yt zdBz1ut Cz1et Δ1z1 ou Az1Δyt zdBz1Δut Cz1et A concepção da representação do processo pela equação 317 é importante para assegurar uma parte integral na lei de controle digital e proporcionar erro nulo em regime entre a saída e a referência para qualquer setpoint e perturbação de carga com comportamento em grau De acordo com a equação 317 diversos tipos de modelos matemáticos lineares discretos mostrados na tabela 33 podem ser utilizados em identificação e controle digital Modelo CARMA Controlled AutoRegressive Moving Average Modelo ARIMA AutoRegressive Integrated Moving Average Modelo CARIMA Controlled AutoRegressive Integrated Moving Average Observações Nas áreas de Identificação e Controle Adaptativo os coeficientes dos polinômios Az¹ Bz¹ e Cz¹ são desconhecidos podem ser constantes ou variantes no tempo e devem ser estimados por técnicas de identificação em uma particular aplicação experimental Malik et al 1991 Shah et al 1991 Levine 1996 Coelho et al 1999 A literatura de Identificação de Sistemas também emprega no lugar dos modelos CAR CARMA CARIMA as notações ARX ARMAX ARIMAX respectivamente onde X denota exogenous ou external input Ljung 1999 36 PROBLEMAS 1 O sistema de controle da figura está sujeito a uma perturbação a partir de uma fonte et Determinar a função de transferência discreta do sistema de acordo com a estrutura CARMA Admitir um período de amostragem de 1 seg e segurador de ordem zero ZOH ZeroOrder Hold Az¹yt Bz¹ut1 Cz¹et 2 Considere o sistema de controle digital ilustrado na figura zerosrootsnumz A função de transferência em malha fechada é dada por Yz z¹1125 0455z¹ 1 0757z¹ 0455z² i Plotar os pólos e zeros no plano complexo z e avaliar a estabilidade em malha fechada ii Qual é a função de transferência em malha aberta Gpz É o processo estável e de fase nãomínima Por quê iii Qual a equação a diferenças do processo polinômios Az¹ Bz¹ e atraso discreto 3 A partir da aproximação numérica da derivada dynTs dt ynTs ynTs Ts Ts obter a correspondente equação a diferenças de segunda ordem para o circuito RLC polinômios Az¹ Bz¹ e atraso de transporte discreto 4 Seja o processo representado pela seguinte função de transferência Gps Ys Us K1Ts s 1Ts1T2s onde K 1 T1 10 seg T2 7 seg e T3 2 seg Considerando que o processo deve ser controlado por computador obter i A representação do processo através da função de transferência discreta supondo o período de amostragem de 1 seg e assumindo o segurador de ordem zero disppólos do sistema discreto ii As alterações na função de transferência z do item anterior se o processo apresentarem atraso de transporte de 4 seg iii A representação do processo CAR em termos da equação a diferenças os polinômios Az¹ e Bz¹ e atraso discreto 5 Para o sistema de controle da figura obter a função de transferência relacionando as transformadas z entre as sequências yt e ut Derivar a equação a diferenças do modelo para Ts 1 seg 6 Considere o projeto dos controladores P PI e PID para manter um nível desejado H2dt no segundo tanque de um sistema de controle de nível de líquido em dois tanques acoplados conforme apresentado na figura 311 O controlador deve comparar o nível atual saída medida H2t yt com o nível desejado H2dt yrt e utilizar a diferença ou erro et H2dt H2t para automaticamente alterar a taxa de fluxo de entrada Q1t ut O correspondente diagrama do sistema está ilustrado na figura 312 polosrootsdenz dnumznumz titleresposta a sequência unitária Controlador ess PO tr tp ts ylabelsaídaxlabelamostras pode ser estabelecido pela variação da tensão ou corrente ou pela abertura ou fechamento de uma válvula de entrada Admitindo uma malha de controle a resposta ao degrau pode ser obtida de acordo com os seguintes passos Mollenkamp 1988 Ogata 2002 i Ajustar o regulador para o modo manual ii Modificar a magnitude da variável de controle acréscimo ou decréscimo iii Registrar plotar a variável de saída do processo Alguns exemplos de resposta ao degrau em malha aberta curva de reação da planta são mostrados na figura 42 A perturbação pode ser definida como Tabela 41 Propriedades do sistema através da resposta ao degrau Resposta Comportamento A estável instável oscilatório atraso fase nãomínima B X X C X X D X X E X X F X Assim a resposta ao degrau é um procedimento adequado para caracterizar de forma preliminar a dinâmica de um processo por meio da simples interpretação gráfica Brosilow e Joseph 2002 Observação Um sistema é conhecido por linear pela determinação da sua resposta ao degrau para diferentes amplitudes no sinal de entrada teorema da superposição um sistema é linear se a forma curva da resposta ao degrau não depende da amplitude de entrada um sinal que tende a afetar adversamente a magnitude da entrada eou e está configurado por três parâmetros ganho estático K constante de tempo τ e atraso de transporte θ Observações A saída de diversos processos químicos exibe respostas ao degrau que podem ser aproximadas pela equação 41 também denominada FOPDT FirstOrder Plus DeadTime O modelo matemático da equação 41 pode ser utilizado na sintonia de diversos controladores PID IMC alocação de pólos Dahlin tempo mínimo entre outros A resposta ao degrau unitário do modelo matemático representado pela equação 41 é ytK1etθ 42 e com base na resposta temporal é possível mensurar pontos para a aplicação de diferentes métodos de estimativa Na literatura existe uma variedade de métodos baseados na resposta do processo ao degrau para identificação de K τ e θ Entre os disponíveis na literatura avaliamse os seguintes métodos apresentados por ZieglerNichols 1942 SundaresanKrishnaswamy 1977 Nishikawa 1984 Smith 1985 e Hägglund 1991 Mollenkamp 1988 Seborg et al 1989 Dorf e Bishop 1995 saída de um sistema Wellstead e Zarrop 1991 Capítulo 4 Métodos clássicos para modelagem de processos 67 A estrutura de um sistema composto do processo e de uma perturbação Capítulo 4 Métodos clássicos para modelagem de processos 67 constante na saída denominada offset está representada na Identificação de sistemas dinâmicos lineares 68 figura 38 isto é conforme ilustra a figura 46 Assim os parâmetros do modelo são dados por K ΔyΔu 412 τ A10368Δy 413 θ t0 τ 414 Figura 46 Método de Nishikawa para a modelagem de processos de primeira ordem Observações Por dependem do traçado da tangente os métodos de ZN e HAG apresentam sensibilidade na presença de ruído imprecisão na sintonia dos parâmetros Os parâmetros estimados são selecionados de tal forma que a soma dos quadrados da diferença entre os valores medidos e os calculados é mínima A função custo no sentido dos mínimos quadrados é dada por J i1N e²i 415 onde N é o número de amostras e ei representa o erro entre a saída real e a saída observada Ljung e Söderström 1983 yt Bz1Az1 utd vt vt constante Com as medidas de resposta ao degrau reais os parâmetros K τ e θ podem variar consideravelmente dependendo das condições de operação do processo isto é da magnitude da mudança do degrau de entrada e da direção da mudança Estas variações podem ser atribuídas às nãolinearidades do processo A seguir o exemplo 41 ilustra a aplicação da técnica de Hägglund na modelagem de um processo de primeira ordem Exemplo 41 De acordo com a curva de reação do processo da figura 47 estimar um modelo paramétrico de primeira ordem pelo método de Hägglund admitir uma entrada degrau unitário yt 0 25 5 75 10 125 15 175 20 225 275 30 tempo min Figura 47 Dinâmica do processo para a aplicação do método de Hägglund A partir da resposta ao degrau unitário identificamse os seguintes parâmetros Ganho valor final 2 K 2 Atraso de transporte θ 12 min Neste tipo de representação da figura 38 é habitual distinguir além Constante de tempo 632 do valor final é 126 126 ocorre em t 435 min τ 435 12 τ 315 min Portanto a função de transferência estimada de primeira ordem é Ȳps 2e12s315s 1 Exemplo 42 Considere o seguinte processo de quinta ordem representado por Gps 2775s5 195s4 141s3 453s2 608s 2775 O teste de resposta ao degrau unitário fornece uma curva de reação conforme ilustra a figura 48 R 0 2 4 6 8 10 tempo seg Figura 48 Resposta ao degrau do sistema da perturbação de carga diferentes tipos de perturbações ou imprecisões que Utilizandose o primeiro método de ZieglerNichols podese identificar um modelo matemático de primeira ordem K 1 θ 06893 e τ 24379 A figura 49 compara as respostas dos processos real modelo de quinta ordem e estimada modelo de primeira ordem para um sinal de entrada degrau unitário Percebese que o modelo identificado não aproxima adequadamente o processo original Ajustes posteriores à aplicação do método de ZieglerNichols são recomendados Obviamente há diferentes formas de extrair informação do processo a partir da curva de reação da resposta ao degrau Alguns exemplos práticos de sistemas de primeira ordem são apresentados na figura 410 afetam a estabilidade do processo controlado tais como perturbação Figura 410 Exemplos de processos de primeira ordem estocástica erro de medição variação paramétrica ou imprecisão na Considerando o servomecanismo de velocidade onde τe τm obtémse a representação em diagrama de blocos dada na figura 411 ou seja figura representação matemática é dada por Vₜs θᵢs Kₘ τₘs 1 onde Kₘ K₁K₂K₃K₄R₂ R₁n Observações Ambos os sistemas apresentados na figura 410 são caracterizados por uma equação diferencial de primeira ordem Um processo de primeira ordem pode ser implementado simulação analógica através de amplificadores operacionais conforme ilustra a figura 412 isto é Figura 412 Implementação analógica da equação 416 estrutura do modelo da planta Quando presente em um sistema de controle e a função de transferência é calculada por E0s leftfracR2R1right2 Eis fracR2R2Cs1 417 o offset é eliminado pela ação integral do controle mas alguns esquemas de A resposta máxima saída no instante do tempo de pico é dada por ytp 1 efractpB 420 e acontece no instante de tempo tp fracpiomeganB 421 controle adaptativo estimam o valor de vt amplitude da perturbação para tr tempo de subida é o intervalo de tempo que o sistema leva para ir de zero a 100 do valor final para sistemas subamortecidos Para sistemas sobreamortecidos utilizase o intervalo de 10 a 90 tp tempo de pico é o tempo gasto para a saída atingir o primeiro pico de overshoot tp fracpiomegansqrt1xi2 422 compensar seu efeito e melhorar a dinâmica de malha fechada da planta i x t2 t1 t3 t1 ii ξ 00805 55470475 x² x 0356 iii f2ξ 07082811ξ se ξ1 f2ξ 26 060 se ξ1 iv ωn f2ξ t3 t1 v f3ξ 0922166ξ vi θ t2 f3ξ ωn vii τ12 ξ ξ² 1 ωn somente se a condição ξ1 é satisfeita A seguir apresentase um exemplo de aplicação do método de Mollenkamp na modelagem de um processo de segunda ordem Exemplo 43 De acordo com a resposta do processo apresentada na figura 414 estimar um modelo paramétrico de segunda ordem pelo método de Mollenkamp variação da entrada degrau de 5 Entretanto as diferentes formas que o sinal vt pode assumir pulso Capítulo 4 Métodos clássicos para modelagem de processos 0 25 50 75 100 125 150 175 200 225 250 275 300 40 41 42 43 44 45 46 47 48 49 50 valor final 495 yt tempo seg Figura 414 Experimentação de um processo para estimação paramétrica 15 de 95 142 t15 50 seg 45 de 95 428 t45 80 seg 75 de 95 712 t75 125 seg Utilizando os instantes t15 t45 e t75 no algoritmo de Mollenkamp obtémse K 19 θ 272 seg τ1 522 seg τ2 199 seg Portanto a função de transferência de segunda ordem estimada é Ğps 19e272s 199s 1522s 1 A tabela 42 ilustra a programação em Matlab do algoritmo de Mollenkamp para o exemplo 43 rampa senoidal ou ruído aleatório O projeto dos algoritmos de Tabela 42 Código em Matlab do exemplo 43 Estimação de um modelo de segunda ordem pelo método de Mollenkamp tlinput valor de t115 t2input valor de t245 t3input valor de t375 xt2t1t3t1 zeta0080555470475x2x0356 if zeta 1 f207082811zeta else f226zeta060 end wmf2t3t1 f30922166zeta zetat2f3wm tall1zetasqrtzeta21wm tall2zetasqrtzeta21wm if zeta 1 dispfator de amortecimento zeta dispfrequência natural wm dispatraso de transporte theta else dispconstante de tempo 1 tall1 dispconstante de tempo 2 tall2 dispatraso de transporte theta end Uma terceira técnica gráfica alternativa para estimar uma função de transferência de segunda ordem é o método de Smith onde sobre a curva de reação do processo são identificados dois instantes de tempo t1 e t2 correspondentes a 20 e 60 do valor final com o aparente atraso de transporte removido Com o auxílio da figura 415 obtémse ξ com o valor de t20t60 e a estimativa de τé calculada a partir de t20t60 x t20t60 controle ótimo depende da natureza da perturbação Na teoria de controle Figura 415 Identificação dos parâmetros tzw t0wr τ ξ pelo método de Smith Para sistemas sobreamortecidos temse Gps Keosτ1s1τ2s1 τ12 τξξ²1 426 enquanto que para sistemas subamortecidos temse Gps Keosτ²s²2ξτs1 427 Um exemplo de sistema de controle de segunda ordem de grande utilização prática em controle de processos está ilustrado na figura 416 estocástico os controladores são desenvolvidos sob a suposição de que as Observações O sistema servomecanismo de posição está caracterizado em malha fechada por uma equação diferencial dominante de segunda ordem A implementação analógica de um processo de segunda ordem através de amplificadores operacionais é apresentada na figura 417 isto é Figura 417 Exemplo de um processo de segunda ordem Admitindo diferentes ganhos 02 03 12 para a implementação do processo de segunda ordem da figura 417 é possível obter comportamentos do tipo sobreamortecido criticamente amortecido e subamortecido Estes resultados estão ilustrados na figura 418 perturbações são estocásticas por natureza A forma importante de 44 MODELAGEM DE PROCESSOS VIA RESPOSTA EM FREQUÊNCIA O método de resposta em frequência na identificação de sistemas lineares está baseado nos diagramas de Bode Neste método entradas do tipo senoidal são aplicadas no sistema e a saída em regime é observada Tanto a razão de magnitude quanto o deslocamento de fase entre a saída e a entrada são avaliados Estas medidas são obtidas em um intervalo de frequências de interesse para o processo particular Se a função de transferência do sistema é Gs então a resposta em frequência é obtida pela substituição de s por jω isto é Gjω Rω jXω 428 Gjω Gωфω Gωejфω 429 Gjω Rω² Xω² 430 фω tg1XωRω 431 Em aplicações práticas devese garantir entradas senoidais em diferentes frequências e medir precisamente nestas frequências as razões de magnitude e o deslocamento de fase O método é aplicável para processos lineares e o procedimento de identificação é offline Além disso o método é aplicável somente para sistemas estáveis desde que a resposta em frequência de sistemas instáveis não pode ser medida na prática Dorf e Bishop 1995 Ogata 2002 Observação Um sistema linear e invariante no tempo estável sujeito a uma entrada senoidal da forma ut Usenωt apresenta em regime uma saída senoidal isto é yt Ysenωt ф Y UGjω 432 com a mesma frequência do sinal de entrada porém com a amplitude e o ângulo de fase da saída diferentes do sinal de entrada perturbação aleatória é aquela associada com mudanças nãoprevisíveis no As diferentes contribuições em termos dos traçados de magnitude e fase dos fatores pólos e zeros comumente encontrados em uma função de transferência são Ganho K 40 MdB 20 20 log Kconstante dB 0 10² 10³ 10⁴ frequência radseg 100 φω 0 100 10² 10³ 10⁴ frequência radseg Pólos ou zeros na origem jωⁿ Admitindo o caso de um pólo na origem N 1 as correspondentes equações de magnitude e fase são 20log ω dB φω 90 90N N 1 em ω 10 20 dB em ω 1 0 dB sistema Estas perturbações são agregadas ao processo e modeladas a partir Capítulo 4 Métodos clássicos para modelagem de processos MdB 0 10¹ 0 50 jω² jω¹ jω⁰ jω¹ jω² 50 10¹ 10⁰ 10¹ 10² 10³ 10⁴ frequência radseg φω 180 90 0 90 180 10¹ 10⁰ 10¹ 10² 10³ 10⁴ Pólos ou zeros sobre o eixo real 1 jωτⁿ Considerando o caso de um pólo simples N 1 então 10log 1 ω²τ² dB φω tg¹ωτ ω 1τ 0 dB ω 1τ 20log ωτ dB ω 1τ 3 dB ω 0 0 ω 90 ω 1τ 45 de uma sequência ruído branco ruído Gaussiano et com média nula e 01τ 1τ 10τ 0 0 φω 100 200 frequência radseg Para o caso de um zero simples N 1 temse o seguinte traçado MdB 100 50 0 01τ 1τ 10τ 100τ 1000τ frequência radseg 100 φω 0 0 100 200 01τ 1τ 10τ 100τ 1000τ frequência radseg Pólos ou zeros complexos conjugados 1 2ξωₙ jω jωωₙ² ⁿ variância σ²e conforme ilustra a figura 39 Admitindo um par de pólos complexos conjugados N 1 então 1 j2ξμ μ²1 μ ωωn 20log Gω 10log 1 μ²² 4ξ²μ² ϕω tg¹ 2ξμ 1 μ² μ 1 magnitude 0 dB ϕω 0 μ 1 magnitude 40 log μ ϕω 180 μ ωωn 1 magnitude na vizinhança de 0 dB depende de ξ ϕω 90 Para ξ 0707 a frequência de ressonância é calculada por ωr ωn 1 2ξ² enquanto que a magnitude máxima de Gjω é Mmax 12ξ1 ξ² 0 ϕω 180 0 01ωωn ωωn 10ωωn frequência radseg Exemplo 44 Considera a função de transferência do processo e o correspondente diagrama de Bode traçado na figura 419 Gps 189 189 s⁴ 59s³ 1203s² 9405s 189 100 Magnitude dB 0 100 200 10¹ 10⁰ 10¹ 10² frequência radseg Figura 419 Diagrama de Bode do sistema de quarta ordem AV12 Pelo diagrama de magnitude identificamse uma frequência de corte em 117 radseg e uma taxa de inclinação de 80 dBdec Isso conduz ao seguinte modelo matemático aproximado ḡps 117⁴s 117⁴ O ganho estático do sistema é unitário e a assintota de alta frequência corta a reta do ganho estático em ω 117 radseg A inclinação de 80 dBdec sugere ao sistema quatro pólos em relação ao número de zeros Para comparar e ter uma ideia da aproximação do modelo estimado perante o sistema real tracemse os diagramas de magnitude e fase conjuntamente na figura 420 45 MODELAGEM DE PROCESSOS VIA RESPOSTA IMPULSIVA Esta seção apresenta um procedimento de estimação dos parâmetros através da resposta impulsiva do sistema sequência de ponderação As condições para a aplicação da técnica são as seguintes Sinha e Kuszta 1983 i A sequência de ponderação está mensurada no universo de interesse da experimentação ii O nível de ruído nas informações do processo deve ser mínimo iii A ordem do modelo do processo deve ser conhecida a priori Considera a função de transferência discreta de malha aberta a ser estimada da planta Gpz de ordem n e descrita por Yz Gpz b0 b1z1 b2z2 bnzn 1 a1z1 a2z2 anzn 433 Admitindo uma sequência impulso unitário na entrada do processo ut1 somente para t0 então a saída é a resposta impulsiva e calculada por ˆHz k0 to hkzk h0 h1z1 h2z2 434 Com as equações 433 e 434 é possível obter a seguinte relação b0 b1z1 b2z2 bnzn h0 h1 a1h0z1 hn Σi1 to n ai hni zn hm Σi1 to n ai hmi zm mn 435 Relacionando os coeficientes de mesma potência da equação 435 obtêmse as n 1 equações com os coeficientes do numerador que podem ser agrupados como b0 1 0 0 0 h0 b1 a1 1 0 0 h1 b2 a2 a1 1 0 h2 hn an an1 an2 1 hn 436 Da equação 436 observase que os coeficientes do polinômio Bz1 são identificados se os coeficientes de Az1 são conhecidos Para determinar os ai i 1 n considere os termos de zn l até z2n Desde que não existem outros termos bi i 1 n no lado esquerdo da equação 435 é possível escrever h1 h2 hn a n hn1 h2 h3 hn1 an1 hn2 hn hn1 h2n1 a1 h2n 437 A solução da equação 437 existe se o rank da matriz contendo a sequência de ponderação é igual à ordem do sistema n Portanto resolvendose primeiramente a equação 437 identificamse os coeficientes do polinômio Az1 pólos de malha aberta para depois com a equação 436 obter os coeficientes do polinômio Bz1 zeros de malha aberta a2 h1 a1 h2 h3 z3 a2 h2 a1 h3 h4 z4 a2 h3 a1 h4 z5 a2 h4 z6 Relacionando os coeficientes de n 1 3 até 2n 4 então a2 h1 h2 h3 a1 h2 h3 h4 Os coeficientes do polinômio Bz1 são obtidos por b0 1 0 0 h0 b1 a1 1 0 h1 b2 a2 a1 1 h2 Exemplo 45 Obter as equações 436 e 437 para um processo caracterizado pela seguinte função de transferência discreta e sequência de ponderação ˆHz b0 b1 z1 b2 z2 h0 h1 z1 h2 z2 h3 z3 h4 z4 1 a1 z1 a2 z2 a3 z3 De acordo com a equação 435 é possível escrever ignorandose os termos acima de 2n b0 b1 z1 b2 z2 h0 h1 a1 h0 z1 a2 h2 a1 h3 z4 𝗚𝘗𝑧𝟕𝟏𝟓𝟕𝟎𝟑𝑙𝑧𝑖𝑛𝑡𝑒𝑟𝑛𝑎𝑡𝑖𝑜𝑛𝑎𝑙𝟔𝟒𝟖𝟕𝟓𝟒𝑡𝑛21𝟐𝟑𝟐𝟓𝟕𝑙𝑧𝟏𝟏𝟕𝟔𝟒𝟎𝟖𝟖𝑡𝑛2𝟎𝟒𝟖𝟕𝟗𝟐𝑖𝑛𝑡𝑒𝑟𝑛𝑎𝑡𝑖𝑜𝑛𝑎𝑙𝑛𝑎𝑛𝑔𝑢𝑎𝑔𝑒 46 PROBLEMAS 1 Aplicar os diferentes métodos de estimação de modelos paramétricos de primeira ordem na identificação dos seguintes processos 𝑮𝒑𝒔𝒆𝟒𝒔𝒔𝟏𝟑 𝑮𝒑𝒔𝟏𝒔𝟏𝟖 2 Considere o sistema de terceira ordem representado por 𝑮𝒑𝑠10𝑠02𝑠105𝑠1 3 Obter a função de transferência estimada para as seguintes curvas de magnitude 4 Considere as medidas da resposta em frequência para uma combinação de servomotor e servoamplificador em um laboratório de controle de processos conforme apresentado a seguir 𝑓𝐻𝑧 20𝑙𝑜𝑔𝑣2𝑣1𝑑𝐵 𝑣2𝑣1𝑔𝑟𝑎𝑢𝑠 01 71 635 02 70 126 03 67 185 04 64 241 05 60 292 06 56 338 07 51 380 08 46 418 09 41 451 10 37 481 12 27 533 15 14 591 20 06 659 25 23 703 30 37 734 40 60 774 50 79 798 70 108 827 100 138 849 5 As amostras da resposta em frequência são apresentadas na tabela Considerando a seguinte função de transferência da planta Gps Kszsps² 2ξωns ωn² plotar as curvas de Bode e identificar os parâmetros K z p ξ ωn do sistema ω radseg 20logG dB G graus 01 032 490 02 110 743 03 204 701 04 293 433 05 372 016 06 496 1077 07 615 3125 10 2581 16851 10 3291 17266 15 3792 17451 20 4498 17625 30 4998 17726 40 5386 17781 50 6590 17890 6 Quando uma entrada degrau foi aplicada em um sistema de segunda ordem a resposta em regime foi medida em 285 A sobreelevação máxima foi de 20 e ocorreu em 25 mseg após a aplicação do degrau Identificar o modelo do sistema 7 As amostras da resposta impulsiva de um processo de segunda ordem são t 0 005 01 015 02 025 03 035 04 ht 0 068 107 127 134 133 128 121 111 Determinar as funções de transferências contínua e discreta da planta 8 Um engenheiro de controle de processos de uma fábrica de motores realiza um ensaio de resposta ao degrau em malha aberta em um servomecanismo de posição Aplicando um degrau na entrada do sistema tensão de alimentação observa a resposta conforme a figura 421 Não identificando os parâmetros do sistema o engenheiro decide implementar um sistema de controle proporcional conforme ilustra a figura 422 realizando sucessivos ensaios de resposta ao degrau em malha fechada com diferentes valores de Kc Após várias respostas decide então avaliar o modelo matemático do sistema como ilustra a figura 422 onde Kc 05 i Qual o modelo matemático obtido pelo engenheiro de processos ii Mostrar que o ganho estático do sistema é unitário e o erro em regime é nulo iii Calcular a resposta do sistema de posição em malha aberta para uma entrada pulso unitário com duração de 10 seg 9 Um engenheiro de controle realiza diferentes ensaios sobre um sistema de controle Na primeira experimentação aplica um degrau unitário Como a resposta é ruidosa o engenheiro somente consegue determinar o valor final da resposta Neste caso y 2 e a derivada na origem de yt é nula Na segunda e na terceira experimentação obtém os seguintes valores em regime u1t sent y1t 225 u2t sen2t y2t 210 i Determinar a função de transferência contínua do sistema a partir das informações dos três ensaios ii Para controlar o sistema implementase uma lei de controle PI digital do tipo ut ut 1 q0et 0905et 1 et yrt yt onde Ts 01 seg e yvt é a sequência de referência do sistema Ajustar o valor de q0 para que o sistema em malha fechada apresente estabililidade Capítulo 5 IDENTIFICAÇÃO DE SISTEMAS REPRESENTADOS POR EQUAÇÕES A DIFERENÇAS Introdução Formalismo histórico dos mínimos quadrados Estimador dos mínimos quadrados nãorecursivo Estimador dos mínimos quadrados recursivo Estimação de processos variantes no tempo Algoritmo de estimação da aproximação estocástica Algoritmo de estimação da variável instrumental Algoritmo de estimação da matriz estendida Problemas 51 INTRODUÇÃO Diferentes métodos para a estimação de parâmetros de modelos lineares discretos são apresentados Ênfase maior é dada ao estimador dos mínimos quadrados uma vez que é a base para o desenvolvimento de outros métodos de identificação Ljung 1999 A qualidade da estimação é dependente da natureza do ruído da estrutura do modelo do tipo de aplicação e da riqueza da informação contida nas medidas Tanto a implementação offline quanto a online do estimador são apresentadas e os aspectos computacionais são avaliados Isermann 1980 Haber e Unbehauen 1990 novo ponto de operação Como em aplicações práticas têmse processos nãolineares o deslocamento do ponto de operação desejado é perigoso do ponto de vista da estabilidade e também pode proporcionar produtos manufaturados fora das especificações de projeto ii o método do teste de resposta ao grau é adequado se o processo sob avaliação contém um baixo nível de ruído variância reduzida o que muitas vezes não acontece em condições práticas industriais iii de forma diferente do teste de resposta ao grau cuja interpretação gráfica é aplicada no conjunto de amostras coletadas da experimentação procedimento para inferir sobre os parâmetros do modelo onde não necessariamente a magnitude da variável do processo medida domina o nível de ruído o estimador dos mínimos quadrados manipula as medidas de entrada e saída nas formas nãoiterativa e iterativa através de algoritmos nãorecursivo e recursivo Considerase a priori que a ordem do modelo é conhecida e que as amostras das medidas de entrada e saída estão disponíveis a cada período de amostragem no universo da experimentação equispaced Os principais objetivos deste capítulo são i familiarizar o leitor com a derivação propriedades e utilização do algoritmo dos mínimos quadrados recursivo MQR quando aplicado na estimação de sistemas lineares discretos SISO SingleInput and SingleOutput com parâmetros desconhecidos eou variantes no tempo ii na inviabilização de aplicação do estimador dos MQR estudar algoritmos alternativos de estimação paramétricos e iii viabilizar a utilização do estimador dos MQR em aplicações de controle adaptativo isto é no contexto de algoritmos de controle autoajustável selftuning control Roffel et al 1989 Åström e Wittenmark 1995 Coelho et al 1999 52 FORMALISMO HISTÓRICO DOS MÍNIMOS QUADRADOS Karl Friedrich Gauss formulou o Princípio dos Mínimos Quadrados ao final do século 18 para prever a trajetória de planetas e cometas a partir das observações realizadas K F Gauss estabeleceu que os parâmetros desconhecidos de um modelo matemático deveriam ser selecionados de modo que o valor mais provável das grandezas desconhecidas é a que minimiza a soma dos quadrados da diferença entre os valores atualmente observados e os valores calculados multiplicados por números que medem o grau de precisão onde quanto mais precisa a medida maior a sua ponderação Ljung e Söderström 1983 53 ESTIMADOR DOS MÍNIMOS QUADRADOS NÃORECURSIVO Considere um processo físico caracterizado por uma entrada ut uma saída yt uma perturbação et e com função de transferência discreta linear da forma Az1yt zdBz1ut et 51 onde Az1 1 a1z1 anzn Bz1 b0 b1z1 bnzn e cuja representação por uma equação a diferenças é yt a1yt 1 a2yt 2 anyt na b0ut d b1ut d 1 bnut d nb et 52 Observações Para o modelo da equação 52 temse na nb 1 parâmetros a estimar Para determinar os ai i 1 na e bj j 0 nb devese utilizar as medidas de entrada e saída do processo O termo et pode representar o erro de modelagem o erro de medição ou o ruído na saída do tipo estocástico determinístico ou offset Definindose o vetor de medidas ϕt de dimensão na nb 1x1 ϕTt yt 1 yt 2 yt na ut d ut d nb 53 e o vetor de parâmetros θt de dimensão na nb 1x1 θTt a1 a2 an b0 b1 bn 54 podese reescrever a equação 52 como yt ϕTtθt et 55 que é denominado modelo de regressão linear Ljung e Söderström 1983 Admitting that N measurements are made sufficient to determine the parameters ai and bj it follows that The weighted least squares estimator also known as the Markov estimator is obtained by minimizing the following criterion Observations J Y φθT WY φθ YTWY YT Wφθ θT φTWY θT φTWφθ W WT Tabela 51 Programação em Matlab do estimador dos mínimos quadrados Identificação de sistemas dinâmicos lineares Capítulo 5 Identificação de sistemas representados por equações a diferenças Na aplicação do estimador dos mínimos quadrados figura 52 todas as medidas devem estar disponíveis a priori para análise e não existe limitação no tempo de processamento do algoritmo identificação offline processo dispositivo de medição computador modelo parâmetros Figura 52 Procedimento de identificação offline A seguir os exemplos 51 e 52 discutem algumas características importantes do estimador dos mínimos quadrados dado pela equação 513 Exemplo 51 Seja um processo sem atraso de transporte representado por yt b₀ut b₁ut 1 et A matriz φ para N amostras é caracterizada por φ ut ut 1 ut 1 ut ut N ut N 1 para os parâmetros estimados b₀ e b₁ A matriz φᵀφ é então calculada por φᵀφ i1 to Nu²i i1 to Nuiui1 i1 to Nui1ui i1 to Nu²i1 Considerações A dimensão de φᵀφ depende do número de parâmetros desconhecidos não do número de amostras Neste caso a matriz é 2x2 Para m parâmetros desconhecidos a matriz tem dimensão mxm A matriz φᵀφ é simétrica de modo que somente a parte triangular superior ou inferior precisa ser calculada Se ut é uma constante ut u₀ então yt b₀ b₁ u₀ et a matriz φᵀφ é dada por φᵀφ Nu₀ 1 1 e assim a matriz φᵀφ é singular e uma solução única mínimos quadrados não pode ser obtida porque φᵀφ não é inversível Para a solução ser obtida então ut deve variar suficientemente para garantir que det φᵀφ 0 quando N cresce Esta condição está geralmente associada com o termo suficientemente excitado No contexto de controle adaptativo é importante que ut mude suficientemente para evitar um rank deficiente para a matriz φᵀφ A seleção de um sinal de entrada particular a ser aplicado na planta pode ser justificada com base na equação da variância assintótica isto é VPM n R N S onde VPM significa a variância nos parâmetros do modelo estimado n é o número de parâmetros do modelo matemático na nb d relativo à equação 52 N é o número de elementos do conjunto de dados da experimentação R é a variância do ruído e S é a variância da entrada A parcela SR razão sinalruído influencia a variabilidade no modelo estimado Condições operacionais altamente ruidosas requerem um grande conjunto de dados eou uma amplitude do sinal de entrada elevada Rivera e Flores 2000 A precisão das estimativas está associada com o tamanho dos elementos da matriz de covariância que por definição é Pt φᵀtφt¹ Variabilidade dos parâmetros estimados Exemplo 52 Obter a função de transferência estimada para um processo representado pelo seguinte modelo matemático yt b₀ut b₁ut 1 utilizando o estimador dos mínimos quadrados procedimento offline Determinar as estimativas para 𝑡 dos parâmetros b₀ e b₁ para N 7 8 14 considerando as medidas de entrada e saída de acordo com a tabela ou seja t 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 ut 1 08 06 04 02 02 04 06 08 1 08 06 04 02 yt 09 25 24 13 10 20 04 14 19 23 24 23 13 12 Para N 7 temse y1 u1 u0 y2 u2 u1 y3 u3 u2 y4 u4 u3 y5 u5 u4 y6 u6 u5 y7 u7 u6 b₀ b₁ A estimativa ótima para b0 e b1 é avaliada por θ ΦᵀΦ¹ΦᵀY e calculada de acordo com Substituindose os valores numéricos da tabela obtémse