Relatório de Modelagem em Espaço de Estado de um Sistema Massa-Mola-Amortecedor

UNIVERSIDADE FEDERAL DE ALAGOAS INSTITUTO DA COMPUTAÇÃO ENGENHARIA DA COMPUTAÇÃO Avaliação AB1 Equipe Anderson Miguel Clemente Santos Ester de Lima Pontes Andrade Guilherme de Oliveira Monteiro Peixoto Lilian Fabrício Marques Neves Professor Ícaro Bezerra Queiroz de Araújo 2 UNIVERSIDADE FEDERAL DE ALAGOAS INSTITUTO DA COMPUTAÇÃO ENGENHARIA DA COMPUTAÇÃO Avaliação AB1 Relatório referente a um sistema físico seu funcionamento modelagem do sistema em espaço de estado descrição da função de transferência de observabilidade e controlabilidade equivalente simulação do modelo e análise como requisito para avaliação da disciplina Sistemas de controle 2 MaceióAL 2022 3 RESUMO Este presente relatório pretende apresentar a implementação de um sistema de terceira ordem realizando a modelagem em Espaço de Estado Sendo assim este projeto modela o sistema massa mola amortecedor que usa duas molas dois carrinhos representativos e dois amortecedores para a construção do sistema 4 SUMÁRIO Introdução 5 Apresentação do Sistema 5 Modelagem em espaço de estado 5 Cálculo da Função de Transferência 8 Simulação do Modelo 9 Análise e comportamento 11 TipoEstabilidade 11 Desempenho do sistema 12 Controlabilidade 13 Observabilidade 13 Implementação de Servo controlador 13 5 1 Introdução Este presente relatório tem como objetivo apresentar a implementação de um sistema de terceira ordem realizando a modelagem em Espaço de Estado Sendo assim este projeto modela o sistema massa mola amortecedor que usa duas molas dois carrinhos representativos e dois amortecedores para a construção do sistema Para a implementação foi utilizado a linguagem Python com auxílio das bibliotecas scipy e control respectivamente para a construção das matrizes e para função de transferência Neste relatório iremos descrever o sistema e seu funcionamento por sequência a modelagem do sistema em espaço de estado a descrição da função de transferência de observabilidade e controlabilidade equivalente apresentaremos a simulação do modelo e por fim traremos a análise e discussão final 2 Apresentação do Sistema Neste projeto iremos utilizar com representação da massa dois carrinhos de rolamento como mostrado na Figura 1 na figura abaixo Figura 1 Representação Esquemática do sistema 3 Modelagem em espaço de estado Considere o sistema apresentado na Figura 1 acima As variáveis de interesse estão na figura e definidas como M1 M2 massas dos carrinhos pt qt posição dos carrinhos ut força externa que está agindo no sistema k1 k2 constantes de mola e b1 b2 coeficientes de amortecimento Os diagramas de corpo livre das massas M1 e M2 são mostrados na figura 2 abaixo Presumimos que o atrito de rolamento dos carrinhos é desprezível 6 Dado o diagrama de corpo livre com forças e direções apropriadas aplicadas usaremos a segunda lei de Newton somatório das forças é igual a massa do objeto multiplicado pela sua aceleração para obter a equação de movimento uma equação para cada massa Para M1 temos Σ𝐹 𝑀1𝑝 𝑢𝑡 𝑘1𝑝 𝑞 𝑏1𝑝 𝑞 𝑀1𝑝 𝑀1𝑝 𝑘1𝑝 𝑏1𝑝 𝑢𝑡 𝑘1𝑞 𝑏1𝑞 𝑒𝑞𝑢𝑎çã𝑜 𝐼 Similarmente para temos 𝑀2 Σ𝐹 𝑀2𝑞 𝑘1𝑝 𝑞 𝑏2𝑝 𝑞 𝑘2𝑞 𝑏2𝑞 𝑀2𝑞 𝑀2𝑞 𝑞𝑏1 𝑏2 𝑞𝑘1 𝑘2 𝑘1𝑝 𝑏1𝑝 𝑒𝑞𝑢𝑎çã𝑜 𝐼𝐼 Onde pt qt aceleração de M1 e M2 respectivamente Agora temos um modelo dado pela equação ordinária de segunda ordem nas equações I e II Podemos começar a desenvolver um modelo de espaços de estados definindo x1t pt x2t qt Denotando as derivadas de x1t e x2t como x3t e x4t respectivamente segue que x3t pt x1t x4t qt x2t Substituindo as variáveis de estado na equação 1 temos 𝑀1𝑝 𝑘𝑝 𝑏𝑝 𝑢𝑡 𝑘1𝑞 𝑏1𝑞 𝑀1𝑥3 𝑘1𝑥1 𝑏1𝑥1 𝑢 𝑘1𝑥2 𝑏1𝑥4 Tirando a derivada de x3 temos 𝑥3 1 𝑀1 𝑢 𝑘1 𝑀1 𝑥2 𝑏1 𝑀1 𝑥4 𝑏1 𝑀1 𝑥3 𝑘1 𝑀1 𝑥1 Similarmente para a equação 2 𝑀2𝑞 𝑞𝑏1 𝑏2 𝑞𝑘1 𝑘2 𝑘1𝑝 𝑏1𝑝 𝑀2𝑥4 𝑥4𝑏1 𝑏2 𝑥2𝑘1 𝑘2 𝑘𝑥1 𝑏1𝑥3 𝑥4 𝑘1 𝑀2 𝑥1 𝑏1 𝑀2 𝑥3 𝑘1𝑘2 𝑀2 𝑥2 𝑏1𝑏2 𝑀2 𝑥4 7 Temos agora a expressão para cada estado e a expressão para cada derivada desses estados podemos proceder para a escrita desses estados em forma matricial essa forma matricial é dado pela expressão 𝑥 𝐴𝑥 𝐵𝑢 Se escolhermos pt como a saída do sistema então 𝑦𝑡 𝐶𝑥 𝐷𝑢 𝑦𝑡 1 0 0 0𝑥 0𝑢 Por fim é necessário o uso do método signalStateSpace que irá receber quatro argumentos que vai representar o sistema dentro do espaço de estado em forma de matriz Abaixo o código utilizado para este fim Variáveis de estados de espaço A nparray 00 00 10 00 00 00 00 10 k1m1 k1m1 b1m1 b1m1 k1m2 k1k2m2 b1m2 b1b2m2 B nparray 00 00 10 00 C nparray10 00 00 00 D 00 sys1 signalStateSpaceA B C D t1 y1 signalstepsys1 pltplott1 y1 8 Para essa modelagem vamos supor que o sistema tenha os seguintes parâmetros m1 m2 1kg k1 k2 1Nm e b1 b2 1Nsm Variáveis físicas massa m1 10 m2 10 coeficiente das molas k1 10 k2 10 coeficiente dos amortecedores b1 10 b2 10 4 Cálculo da Função de Transferência Utilizando as matrizes modeladas em Espaço de Estados podemos encontrar a função de transferência equivalente utilizando o método signalss2tfABCD da biblioteca scipy tfsys1 signalss2tfABCD convertendo Espaço de Estados para função de transferencia printtfsys10 indices dos zeros printtfsys11 indices dos polos output 000000000e00 399680289e15 100000000e00 200000000e00 200000000e00 1 3 4 2 1 Como a saída não tem uma formatação amigável podemos utilizar o método ctrltf da biblioteca control para exibir a função de transferência de forma mais apropriada Função de transferencia visualizavel system ctrltfnparray 0 0 1 2 2 nparray 1 3 4 2 1 printsystem s2 2 s 2 9 s4 3 s3 4 s2 2 s 1 5 Simulação do Modelo Para a simulação do modelo os sinais de entrada foram definidos da seguinte forma 𝐷𝑒𝑔𝑟𝑎𝑢 𝑢 𝑘 𝑡 0 0 𝑡 0 𝑅𝑎𝑚𝑝𝑎 𝑟𝑡 𝑡 𝑃𝑎𝑟á𝑏𝑜𝑙𝑎 𝑝𝑡 𝑡 2 Também foi definido um vetor de tempo com período de 60 segundos time nparange0 60 001 tempo ustep npfulllentime 1 entrada do degrau uramp npfulllentime time entrada da rampa uparabola npfulllentime time2 entrada da parabola printustep printuramp printuparabola 1 1 1 1 1 1 0000e00 1000e02 2000e02 5997e01 5998e01 5999e01 00000000e00 10000000e04 40000000e04 35964009e03 35976004e03 35988001e03 Simulamos a saída do sistema utilizando o método signallsim do scipy que recebe como parâmetros o sistema dado em variáveis de estado a entrada calculada anteriormente e o tempo também definido anteriormente tout y1 x1 signallsimsys1 ustep time resposta ao degral tout y2 x2 signallsimsys1 uramp time resposta a rampa tout y3 x3 signallsimsys1 uparabola time resposta a parabola Definindo uma função para plotar os resultados def outPutGrapht y signal yLabel pltplott signal orange labelSinal pltplott y blue labelyLabel pltgridalpha03 pltxlabelt pltylabelAmplitude pltlegend 10 pltshow print outPutGraphtime y1 ustep Degrau outPutGraphtime y2 uramp Rampa outPutGraphtime y3 uparabola Parabola 11 6 Análise e comportamento 61 TipoEstabilidade Podemos observar o comportamento do sistema verificando os autovalores da matriz A que representam os pólos desse sistema poles nplinalgeigA0 printpoles 130901699095105652j 130901699095105652j 019098301058778525j 019098301058778525j fsf signalplacepolesA B poles methodYT t nplinspace0 2nppi 401 pltplotnpcost npsint k unit circle pltplotfsfrequestedpolesreal fsfrequestedpolesimag wo labelDesired pltplotfsfcomputedpolesreal fsfcomputedpolesimag bxlabelPlaced pltgrid pltaxisimage pltaxis15 15 11 11 pltlegendbboxtoanchor105 1 loc2 numpoints1 12 Um sistema é Tipo 0 pois não apresenta pólos na origem Portanto calculando todos as constantes de erro estático 62 Desempenho do sistema ctrlstepinfosystem RiseTime 22144577273071495 SettlingTime 18137463290325226 SettlingMin 1754119178357359 SettlingMax 2682536336352232 Overshoot 34126816817611605 Undershoot 0 Peak 2682536336352232 PeakTime 53779687663173625 SteadyStateValue 20 13 7 Controlabilidade MControlabilidade ctrlmatlabctrbAB printMControlabilidade 0 1 1 1 0 0 1 2 1 1 1 1 0 1 2 2 postoControl nplinalgmatrixrankMControlabilidade printpostoControl 0 1 1 1 0 0 1 2 1 1 1 1 0 1 2 2 Uma vez que o posto da matriz de controlabilidade é igual a ordem do sistema podemos concluir que o sistema é completamente controlável 8 Observabilidade obs ctrlmatlabobsvAC printobs 1 0 0 0 0 0 1 0 1 1 1 1 2 3 1 2 postoObs nplinalgmatrixrankobs printpostoObs 4 Uma vez que o posto da matriz de observabilidade é igual a ordem do sistema podemos concluir que o sistema é completamente observável 9 Implementação de Servo controlador Dado a representação do sistema em espaço de estados podemos modificar seu comportamento através de um controlador que atua sobre a matriz de estados A de modo que os autovalores e polos são realocados para uma posição préestabelecida Como o sistema deve ser controlado para uma posição fixa utilizamos o servocontrolador do tipo 0 14 Podemos definir a partir do diagrama de blocos do controlador que a entrada do sistema pode ser definida como eq I 𝑢 𝑟𝑘𝑝 𝑘𝑥 Podemos inserir a equação I na equação das variáveis de estados então 𝑥 𝐴𝑥 𝐵𝑢 𝑥 𝐴𝑥 𝐵𝑟𝑘𝑝 𝑘𝑥 𝑥 𝐴𝑥 𝐵𝑟𝑘𝑝 𝐵𝑘𝑥 𝑥 𝐴 𝐵𝑘𝑥 𝐵𝑟𝑘𝑝 Com base nesta equação a implementação do controlador está definida nos passos 1 Definir a entrada do sistema e o tempo 2 Escolher novos polos para o cálculo da matriz de ganho 3 Encontrar a matriz A em malha fechada 4 Encontrar o Kp 5 Calcular a matriz B para o sistema em malha fechada Definimos o sistema em malha aberta sysopen matlabssA B C D Atribuímos os vetores de tempo e entrada do sistema O tempo está definido de 0 a 150 e o vetor de entrada é uma concatenação de 4 sinais constantes definidos em 5 10 15 e 20 ou seja o sinal de entrada escolhido foi um degrau que varia em 4 valores diferentes com o passar do tempo tempo nparange0 150 001 u1 npfullintlentempo4 5 u2 npfullintlentempo4 10 u3 npfullintlentempo4 15 u4 npfullintlentempo4 20 u npconcatenateu1 u2 u3 u4 Em seguida partindo do vetor de tempo e do vetor que representa o sinal de entrada obtemos a resposta do sistema Calculamos a saída do sistema 15 yo to xo matlablsimsysopen u tempo Essa saída representa o sistema em malha aberta com as entradas degraus em 5 10 15 e 20 Escolhendo arbitrariamente as posições dos pólos desejados 51j e 101j Os polos escolhidos para o sistema em malha fechada obtidos com a análise da frequência natural amortecida e não amortecida constante de atenuação e coeficiente de amortecimento apresentam um comportamento com menor sobressalto P 51j51j101j101j A matriz de ganho K pode ser calculada através da Fórmula de Ackerman definida por onde e é um vetor linha 1xn é a matriz de controlabilidade do 𝑘 𝑒𝑀𝑐 1𝑃𝐴 𝑀𝑐 1 sistema e PA é uma matriz formada a partir do polinômio característico desejado para a malha fechada Neste projeto foi utilizado o método placeA B P da biblioteca controlmatlab que utiliza o algoritmo de Tits and Yang para calcular a matriz de ganho K matlabplaceA B P printK 817 892 9 1684 Com a matriz de ganho podemos calcular a matriz A em malha fechada dado que 𝐴𝑐𝑙 𝐴 𝐵𝑘 Acl A npdotB K printAcl 0 0 1 0 0 0 0 1 917 792 10 1584 1 2 1 2 Polos do sistema em malha fechada 16 printnplinalgeigAcl0 51j 51j 101j 101j Recalculando a saída do sistema em malha fechada com o novo valor da matriz de espaço de estados sysclose matlabssAclBCD yc tc xc matlablsimsysclose u tempo As posições escolhidas para os polos mudaram a dinâmica do sistema Resta garantir que o erro em regime permanente seja o menor possível Calculando o coeficiente de erro 𝐾𝑝 Kdc matlabdcgainsysclose Kp 1Kdc printKp 13129999999999988 Após calcularmos o ganho do sistema podemos calcular o escalar que multiplica a entrada para compensar o resultado do sistema na saída e atribuir esse valor multiplicativo ao módulo do sistema 𝐾𝑝 Calculando a matriz B em malha fechada Bcl npdotKp B Finalmente o sistema pode ser remontado com todas as variáveis encontradas anteriormente sendo os novos parâmetros Acl para as variáveis de estados e Bcl a nova variável de entrada do sistema sysclosekp matlabssAcl Bcl C D Resposta do sistema em malha fechada 17 Comparando as saídas em malha aberta e malha fechada ctrlstepinfosysopen RiseTime 2214457727307151 SettlingTime 18137463290325236 SettlingMin 17541191783573615 SettlingMax 26825363363522428 Overshoot 3412681681761214 Undershoot 0 Peak 26825363363522428 PeakTime 5377968766317366 SteadyStateValue 20 ctrlstepinfosysclosekp RiseTime 3219949671876136 SettlingTime 5366582786460226 SettlingMin 09017957808891607 SettlingMax 09999999999999998 Overshoot 0 Undershoot 0 Peak 09948534242170397 PeakTime 6907755278982137 SteadyStateValue 09999999999999998 18 Referências 1 What is Pole PlacementSite YouTube Online Acessado em 09102022 httpswwwyoutubecomwatchvFXSpHy8LvmY 2 Controle por realimentação dos estados reguladores de estado Site Universidade de São Paulo Online Acessado em 09102022 httpsitespoliuspbrdpmr24001220ReguladorControle20por20realimentaC3A7C3A 3o20dos20estadospdf 4 Place Site Universidade do MathWorks Online Acessado em 15102022 httpsitespoliuspbrdpmr24001220ReguladorControle20por20realimentaC3A7C3A 3o20dos20estadospdf 5 RICHARD C DORF BISHOP Robert H Modern control systems Pearson Prentice Hall 2008