Introdução sobre Identificação Utilizando Modelos de Hammerstein

Ainda que a maioria das contribuições para o projeto de controladores sejam baseadas em modelos lineares, tipicamente processos industriais são inerentemente não lineares, tais como controle de nível em tanques cônicos, controle de temperatura por impedâncias térmicas e controle de vazão por meio de válvulas. Uma metodologia de linearização aplicada nesses casos é a construção de um modelo para cada ponto de operação. Contudo, quando o sistema é fortemente não linear e o ponto de operação muda ao longo de uma ampla região, é difícil representar adequadamente um dado sistema por meio dessa abordagem, dada a grande quantidade de modelos necessária.

Para obter uma boa representação dos sistemas descritos, é necessário modelá-los usando técnicas não lineares, visto que as técnicas lineares fazem aproximações que resultam em modelos pouco confiáveis. Uma classe de modelos não lineares frequentemente estudados, são os chamados modelos não lineares de blocos interconectados. Esses modelos consistem da interconexão de um sistema LIT (Linear Invariante no Tempo) e uma não linearidade estática. Dentro dessas classes, as duas estruturas mais comuns são a de Wiener, em que um bloco dinâmico linear precede um bloco estático não linear. A outra é a estrutura de Hammerstein apresentada na Figura (1), que possui os mesmos blocos, porém com a ordem invertida. O sinal de entrada é representado pela letra u e o de saída por y. A não linearidade estática mapeia a entrada u para a variável intermediária v. E o bloco dinâmico linear mapeia a sequência intermediária v para a sequência de saída do modelo y. O símbolo chapéu (^) indica que os sinais são estimativas. w e v são ruídos de processo e medição respectivamente.

Figura 1: Modelo de Hammerstein em malha aberta.

Ainda possui uma terceira, que combina essas duas abordagens utilizando um bloco dinâmico linear entre dois estáticos não lineares. Modelos de Wiener produzem resultados mais acurados quando o sistema tem dinâmica variável com o ponto de operação, ao passo que o modelo de Hammerstein se adequa melhor quando o sistema tem ganhos variáveis. Essas estruturas foram usadas com sucesso para representar sistemas não lineares em diversas aplicações práticas na área de processos químicos, biológicos e de controle.

O modelo de Hammerstein é simples de ser compreendido e obtido. Por isso será utilizado para exemplificar a utilidade dos modelos de blocos interconectados no contexto de identificação de um sistema dinâmico não linear. 

Com o objetivo de fazer identificação, o primeiro passo consiste em coletar dados da planta que se deseja identificar. Caso seja possível fazer um teste projetado para coletar os dados de forma a extrair informações de interesse, ótimo! Caso não seja, deve-se trabalhar com os dados que tiver disponível. No exemplo dado a seguir os sinais que deveriam ser coletados da planta numa situação prática serão gerados por um modelo, ou seja, será possível projetar um sinal de teste capaz de permitir a extração de informação da "planta".

Exemplo didático: O exemplo apresentará uma maneira de identificar um modelo de Hammerstein. De forma a simplificar, a planta será descrita matematicamente por uma não linearidade estática polinomial, mais especificamente por:

Veja que não tem o chapéu (^) no v visto que se trata de um sinal da planta. Mas perceba que numa planta real esse sinal não existe, apenas o sinal de entrada e o de saída. A dinâmica da planta será descrita pelo modelo ARX de primeira ordem:
Definida a estrutura da planta, os sinais de entrada devem ser projetados. A dinâmica é a mesma independente do ponto de operação, porém o ganho estático não linear muda com o ponto de operação. Por exemplo, quando a entrada é 2 o sinal intermediário vale 4, porém quando a entrada é 10 o sinal intermediário vale 100. Dessa forma o teste da dinâmica pode ser feito em uma região de operação pequena enquanto o testes estático deve ser feito em uma região ampla. 

Quanto deve ser a amplitude do teste estático? Depende da aplicação. Para o exemplo em questão, vamos supor que a planta só permite entrada variando de 0 a 10. Então trabalharemos no limite para o teste estático. 
Como deve ser o sinal do teste estático? Dada a intenção de extrair a característica estática da planta vamos aplicar um sinal que permita mapear a relação entre u e v.  Para facilitar a visualização o sinal de entrada será variado de forma gradual em degraus subsequentes cada vez de amplitudes maiores. A granularidade dos degraus depende da curva estática. No exemplo atual utilizaremos degraus variando de 0 a 10 de 1 em 1 unidade.
Quanto deve ser a amplitude do teste dinâmico? Como dito, é comum que a amplitude do teste dinâmico seja menor, pois causa menos variações na planta (num caso prático) e trabalha em torno de uma faixa de operação em que o sistema responde aproximadamente linear. Dessa forma a amplitude ficará do sinal de teste será limitada entre 0 e 1.
Como deve ser o sinal do teste dinâmico? O sinal deve ser tal que exprima as características dinâmicas da planta na saída. Pensando em termos de frequência, matematicamente, a planta pode ser representada em uma faixa de frequências e se essa faixa for excitada a saída representará bem a dinâmica da planta. Mas qual faixa de frequências? Sem ter que saber essa resposta pode-se inserir na entrada um sinal que excite todas frequências. Um sinal que atende esse requisito é o ruído branco dada sua característica aleatória. Do mesmo modo, é comum empregar o PRBS que é um sinal com amplitudes variando entre dois níveis cuja mudança de nível é aleatória. Será empregado um sinal PRBS variando de 0 a 1. Para mais detalhes do sinal PRBS ver PRBS MATLAB.

Teste estático: Aplicar o sinal de teste estático e coletar a saída em MATLAB: 
deg = 0:1:10;
tempo_degrau = 30;
N = 11*tempo_degrau;
u = reshape(repmat(deg,tempo_degrau,1),1,N);
v = u.^2;
y = zeros(N,1);
for k = 2:N
    y(k) = 0.8*y(k-1) + 0.3*v(k-1);
end

Figura 2: Teste estático.


Em primeiro lugar, perceba que numa situação prática só se tem acesso à entrada e a saída, portanto o ganho que aparece no gráfico é o ganho total da planta (ganho DC da dinâmica junto com o ganho estático não linear). Dessa forma o que se pode fazer é assumir que o modelo dinâmico a ser identificado terá ganho DC unitário e portanto o sinal intermediário em regime permanente será tomado como o próprio sinal de saída. 
Em segundo lugar, perceba que o "tempo de degrau" foi selecionado como 30 iterações. Ele deve ser tal que o sistema possa responder e estabilizar, como apresenta-se na Figura 2, uma vez que o que interessa no teste são os dados em regime permanente. 
Dito isso coletando apenas os sinais em regime permanente de interesse, ou seja as amostras finais de cada degrau, obtém-se da Figura 2 o seguinte conjunto de pontos que caracteriza a curva estática não linear da planta: 
          u = [0   1    2     3      4    5      6      7     8      9        10]    
y = v = [0  1,5  6  13,5  24  37,5  54  73,5  96  121,5  148]

Pode-se aproximar uma curva não linear para esse conjunto de dados utilizando o comando polyfit do MATLAB:
u = [0   1   2    3   4    5    6    7    8     9     10];   
v = [0  1.5  6  13.5  24  37.5  54  73.5  96  121.5  148];
poli = polyfit(u,v,2);
que resulta no seguinte polinômio:
Ou seja, foi identificado um polinômio de segunda ordem para representar a não linearidade estática. A escolha do termo de segunda ordem é baseada na aparência da curva da Figura 2 (observando apenas os dados em regime permanente).

Teste dinâmico: aplicar o sinal PRBS e coletar a saída em MATLAB:
Nd = 256;
u = 4*prbs(Nd,10,1);
v = u.^2;
y = zeros(Nd,1);
for k = 2:Nd
    y(k) = 0.8*y(k-1) + 0.3*v(k-1);
end
O Nd indica número de dados disponível para identificação que deseja-se gerar. O sinal PRBS é gerado pela função prbs() e em seguida inserido na função estática não linear que por sua vez é entrada do modelo dinâmico simulado no laço for. O sinal PRBS e a saída são apresentados na Figura 3.
Figura 3: Sinal de teste e saída coletada correspondente do ensaio dinâmico.


Utilizando o algoritmo de mínimos quadrados 
Psi = [y(1:end-1) u(1:end-1)']; %construção do regressor
theta = inv(Psi'*Psi)*Psi'*y(2:end); %estimação dos parâmetros
Resultando no modelo ARX:

que completa a identificação do modelo de Hammerstein.  Então perceba que para o sistema em questão o ganho não linear estático está sempre com um ganho de 1,5 vezes o que deveria na saída do sistema, o motivo é que o ganho DC do sistema original é exatamente 1,5. Pra verificar isso observe esses sinais que já foram mostrados e será repetido aqui por conveniência:
           u = [0   1    2     3      4    5      6      7     8      9        10]    
y = v = [0  1,5  6  13,5  24  37,5  54  73,5  96  121,5  148]
Mas pensando que o a lei não linear é quadrática então 4^2 deveria resultar em 16, porém o valor coletado do ensaio foi 24. O motivo é que o 16 entra no sistema dinâmico e o ganho DC da parte dinâmica atua sobre ele, multiplicando-o por um fator de 1,5. O ganho DC do sistema original do exemplo pode ser encontrado ao fazer a transformada z da função e fazer o limite de z tendendo a 1. 
O mesmo vale para o modelo identificado. Porém seu ganho DC será encontrado como 0,6 e o restante compensado pelo sinal intermediário como será visto na Figura 4.

O próximo passo é validar o modelo identificado por simulação livre, para isso deve-se gerar um novo conjunto de dados tal como o que foi apresentado na Figura 3, utilizando um sinal PRBS diferente. O teste será realizado nas mesmas condições anteriores, ou seja a amplitude do sinal PRBS será entre 0 e 4 e se trata de uma realização diferente do sinal de entrada. Isso é necessário para verificar a capacidade de generalização do modelo.
u = 4*prbs(Nd,10,1);
v_hat = polyval(poli,u)*(1-theta(1))/theta(2);
y_hat = zeros(Nd,1);
for k = 2:Nd
    y_hat(k) = theta(1)*y_hat(k-1) + theta(2)*v_hat(k-1);
end

v = u.^2;
y = zeros(Nd,1);
for k = 2:Nd
    y(k) = 0.8*y(k-1) + 0.3*v(k-1);
end

Figura 4: A esquerda apresenta-se o sinal intermediário e a direita o sinal de saída.

Pode-se notar que o sinal intermediário é diferente do original em termos de amplitude, porém ela é compensada pelo ganho que o modelo dinâmico insere no sistema identificado de modo que o sinal de saída acompanha o sinal do sistema original, ou seja, comparando entrada/saída o modelo de Hammerstein identificado é capaz de descrever o sistema do exemplo.


Autor: Luís Henrique dos Santos

SANTOS, L. H. Introdução sobre identificação utilizando modelos de Hammerstein. Blog: Identificação e Controle de Sistemas Lineares e Não Lineares. ago. de 2022. Disponível em: https://identificacao-e-controle.blogspot.com/2022/07/modelos-de-blocos-interconectados.html

Comentários

Postagens mais visitadas deste blog

Aplicativo de Identificação de Sistemas Utilizando MATLAB App Designer