PID digital – Método de cálculo numérico

Introdução

No artigo técnico Controlador PID digital: Uma modelagem prática para microcontroladores – Parte 1, é  apresentado com detalhes o controlador PID, sua arquitetura típica, o desenvolvimento das equações para o controle e alguns exemplos de como se pode desenvolver um controle PID. Nesse artigo serão retomados alguns conceitos sobre esse tipo de controle e serão apresentadas formas de cálculo numérico para implementar a parte diferencial e a integral do controlador. Essa forma alternativa de cálculo casa muito bem com a aplicação em microcontroladores ou processadores de 8 bits ou outros relativamente lentos. Os resultados que serão apresentados a seguir, foram largamente aplicados nos primeiros controladores PID digitais e certamente ainda são úteis nos dias de hoje, por sua simplicidade e velocidade de processamento.

Arquitetura típica de um controlador PID

Um sistema de controle típico de PID em malha fechada tem o aspecto mostrado na Figura 1. Nessa figura pode-se observar uma entrada, referente a um sinal de referência ou setpoint, um bloco onde essa referência é subtraída de um sinal proveniente do processo, em geral medido de alguma forma por um sensor, transdutor etc. Essa subtração produz um sinal de erro que entra no controlador PID propriamente dito, e a saída do controlador, que atua sobre o processo que está sendo controlado, é aplicada para conduzir o processo ao estado desejado.

.

.

Na prática, o controlador PID digital destacado na Figura 1 em vermelho, pode ser representado como mostrado na Figura 2. O erro é amostrado, produzindo assim um conjunto de amostras discretizadas do sinal de erro, que é conduzido a um conversor analógico / digital (A/D), depois são calculados os valores da derivada, da integral e do valor proporcional, em seguida esses valores calculados são somados e conduzidos a um conversor digital / analógico (D/A).  

O desenho da Figura 2 representa o fluxo do sinal de erro para o controlador PID de forma a simplificar o desenvolvimento dos cálculos que serão apresentados a seguir. É evidente que num controlador real, os sinais de SP (setpoint) e PV (present value) também são digitalizados e subtraídos para que se obtenha o erro. O importante é fixar o conceito de que o erro deve ser entregue ao controlador numa frequência ou taxa fixa. Isso é requisito para o desenvolvimento que será realizado a seguir. Para que o sistema funcione bem, também é necessário que a atualização do sinal de saída pelo conversor D/A seja feito na mesma frequência do amostrador de entrada. Na Figura 3 é mostrado um sinal analógico como exemplo de sinal de erro e são destacadas as amostras realizadas num período fixo (Δ).

.

.

A seguir será desenvolvido um pouco de teoria, o suficiente para que se possa saber de onde foram tiradas as equações que resultarão desse desenvolvimento. Se você preferir, pode consultar direto os quadros resumo que apresentam as equações resultantes desses cálculos.

 .

Séries de Taylor

Brook Taylor elaborou uma teoria a respeito de funções matemáticas em 1715, conhecida como teorema de Taylor ou séries de Taylor. Esse teorema estabelece que: sendo f  uma função  infinitamente diferenciável, ela pode ser representada pela seguinte fórmula:

(1)

onde:

• 

    representa a k-ésima derivada da função avaliada no ponto x₀;

• 

            é um ponto arbitrário de referência;

• 

            é o fatorial de k.

Se a somatória acima convergir, a função pode ser representada por uma série finita de termos. Como por exemplo a expansão da série (1): (2)

Obs: Quando x₀ = 0, a série (2) é conhecida como série de McLaurin.

A seguir serão utilizadas as séries de Taylor para deduzir formas simples para o cálculo da primeira derivada de uma função e a integral dessa função.

Derivação Numérica

Seja definida uma função f(x), representada por uma série de Taylor, tal como mostrada em (1), onde x = i é uma amostra arbitrária de índice i tal como mostrado na Figura 3. O valor de f(x) será batizado de f_i. Seja f_(i-1) o valor dessa mesma função calculado no índice (i-1). A distância entre essas amostras é de Δ, correspondente a uma unidade de tempo relacionada com a frequência ou taxa de amostragem. Daí obtemos:

(3)

Para avaliar a função em (i-1), vamos substituir i por (i-Δ). Realizando a substituição e rearranjando a equação obtida, chega-se na equação (4).

(4)

Escrevendo essa somatória na forma expandida, obtém-se:

 (5)

Se desprezarmos os termos de segunda ordem e superiores e rearranjando, a equação fica assim:

 (6),  e o erro =

.

A expressão (6) é conhecida como a derivada com 2 pontos.

.

É evidente que nesse cálculo, o erro é relativamente grande. Pode-se estender o raciocínio desenvolvido acima para um cálculo mais preciso, que pode ser obtido com 3, 4 ou mais amostras do sinal.

.

Se substituirmos i por (i- 2*Δ) e for desenvolvido o mesmo raciocínio que foi utilizado para i-1 poderemos levantar a equação de f_(i-2) em função de f_i.

  (7)

Usando a equação (5) e rearranjando, obtém-se:

.

 (8), onde o erro = 

.

Essa é a equação de cálculo de derivada com 3 pontos. De maneira análoga pode-se desenvolver a derivada de 4 pontos, mostrada na equação (9).

.

(9), onde o erro = 

.

A seguir será apresentado um quadro resumo com as fórmulas desenvolvidas até aqui.

.

Quadro resumo da diferenciação digital

Aqui será adotada uma nova notação, mais frequentemente utilizada em processamento digital de sinais (DSP) e em controle digital. Isso facilitará o entendimento de quem está habituado a essa notação e introduz a notação para quem ainda não a conhece. A notação é assim:

• y = saída (nesse caso a derivada da função);
• x = amostra atual da função;
• z^-n = n-ésima amostra anterior;
• Δ = intervalo de tempo entre amostras.

Resumo das fórmulas referentes à diferenciação digital:

.

=====================================================================

• Derivada calculada com 2 pontos:

(10)

 —————————————————————————————————————————

• Derivada calculada com 3 pontos:

(11)

————————————————————————————————————————–

• Derivada calculada com 4 pontos:

(12)

===================================================================

.

Os resultados obtidos nas equações acima são facilmente implementáveis em microprocessadores.

.

Integração Numérica

Para desenvolver as equações referentes à integração numérica, o raciocínio é o mesmo já desenvolvido para a diferenciação. Voltando à Figura 3, onde temos representada uma forma de onda arbitrária amostrada em intervalos regulares de Δ segundos, pode-se definir que:

(13)

Seja o segundo termo da equação (13):

(14)

Desenvolvendo a série para I_i-1 se obtém uma equação semelhante à equação (5):

(15)

Rearranjando, obtém-se:

(16)

Desenvolvendo:

(17)

 Substituindo:

(18)

desprezando os termos de segunda ordem ou superiores chega-se à equação para o cálculo da integral com 2 pontos.

.

(19), onde o erro = 

.

A fórmula (19) é conhecida como “regra do trapézio”.

.

O cálculo da integral pode ser estendido para mais pontos, para aumentar a precisão do resultado. A seguir, a equação para o cálculo com 3 pontos (20) e 4 pontos (21).

.

 (20), onde o erro = 

.

(21), onde o erro = 

.

Quadro resumo da integração digital

.

Aqui também será adotada a nova notação, mais frequentemente utilizada em processamento digital de sinais (DSP) e em controle digital.  A notação é assim:

• y = saída (nesse caso a integral da função);
• x = amostra atual da função;
• z^-n = n-ésima amostra anterior;
• Δ = intervalo de tempo entre amostras.

Resumo das fórmulas referentes à integração digital:

=====================================================================

• Integral calculada com 2 pontos:

(22)

—————————————————————————————————————————

• Integral calculada com 3 pontos:

(23)

————————————————————————————————————————–

• Integral calculada com 4 pontos:

(24)

===================================================================

Observe que as fórmulas acima calculam a integral incremental entre uma amostra e a anterior. Para que seja realizado um cálculo de integral de fato é necessário acumular os resultados.

Obtidos os resultados  acima, é interessante verificar se eles realmente funcionam. Para isso eles serão simulados num programa para cálculos matemáticos conhecido por Scilab. Serão mostrados os respectivos códigos e os resultados. Para que você possa testar e explorar na sua própria máquina será necessário instalar o Scilab nela. Baixe o programa aqui: SCILAB vers. 5.5.0. O programa é gratuito e tem versões para Windows, GNU/Linux e Mac OS X.

Inicie o programa e digite o código a seguir. Se você preferir, você pode marcar o código, copiá-lo (Ctrl-C) e colá-lo no Scilab (Ctrl-V). Para iniciar as operações basta teclar o “Enter”.

.

Cálculo de Integrais numéricas utilizando o Scilab

.

// Inicializar as variáveis e a forma de onda para os testes

nCiclos = 10; // Define o núemro de ciclos da forma de onda
nAmostras = 40; // Define o número de amostras por ciclo

// Gera o vetor de referência para a produção da forma de onda

t = (1:2*%pi/nAmostras:nCiclos * 2 * %pi + (1-2*%pi/nAmostras)); 
nvErro = squarewave(t); // Gera uma onda quadrada com 400 amostras

// Gera o gráfico para essa forma de onda

a0 = figure(0);
a0.background = -2;

plot(nvErro);

// Gera o grid e as legendas

xtitle('Forma de onda do Erro para testes', 'Indice das amostras', 'Amplitude');
set(gca(), "data_bounds", [0,400,-1.2,1.2]);
set(gca(), "grid", [1,1]);

Os comandos acima geram um gráfico parecido com o apresentado na Figura 4.

.

Como não há formas de onda triangulares prontas no Scilab, será utilizada a forma de onda quadrada, calculada sua integral para formar a onda triangular e calculada a diferencial da onda triangular para formar a onda quadrada novamente. Digite o seguinte código no Scilab (ou copie e cole):

.

// Inicializa as variáveis para o cálculo da integral com 2 pontos

Int2p = zeros(1, size(nvErro,2));
Int2pAcum = 0;
Z1 = 0;
Delta = 1;

// Calcula a Integral com 2 pontos


for nIndice = 1:size(nvErro,2)
    Int2p(nIndice) = ((nvErro(nIndice) + Z1)* Delta/2) + Int2pAcum;
    Int2pAcum = Int2p(nIndice);
    Z1 = nvErro(nIndice);
end;

// Apresenta a integral no gráfico

a1 = figure(1);
a1.background = -2;

plot(Int2p, 'm');

// Gera o grid e as legendas

set(gca(), "background", [-2]); 
xtitle('Integral da forma de onda calculada com 2 pontos', 'Indice das amostras', 'Amplitude'); 
set(gca(), "data_bounds", [0,400,-10,20]); 
set(gca(), "grid", [1,1]);

//==========================================================================

// Inicializa as variáveis para o cálculo da integral com 3 pontos

Int3p = zeros(1, size(nvErro,2));
Int3pAcum = 0;
Z1 = 0;
Z2 = 0;
Delta = 1;

// Calcula a Integral com 3 pontos

for nIndice = 1:size(nvErro,2)
     Int3p(nIndice) = ((5*nvErro(nIndice) + 8*Z1- Z2)* Delta/12) + Int3pAcum;
     Int3pAcum = Int3p(nIndice);
     Z2 = Z1;
     Z1 = nvErro(nIndice);
end;

// Gera gráfico

a2=figure(2);
a2.background = -2;

plot(Int3p, 'b');

// Gera o grid e as legendas

set(gca(), "background", [-2]); 
xtitle('Integral da forma de onda calculada com 3 pontos', 'Indice das amostras', 'Amplitude'); 
set(gca(), "data_bounds", [0,400,-10,20]); 
set(gca(), "grid", [1,1]);

//========================================================================================================= // Inicializa as variáveis para o cálculo da integral com 4 ponto

Int4p = zeros(1, size(nvErro,2));
Int4pAcum = 0; 
Z1 = 0; 
Z2 = 0; 
Z3 = 0; 
Delta = 1; 

// Calcula a Integral com 4 pontos 

for nIndice = 1:size(nvErro,2) 
    Int4p(nIndice) = ((9*nvErro(nIndice) + 19*Z1- 5*Z2 + Z3)* Delta/24) + Int4pAcum;
    Int4pAcum = Int4p(nIndice); 
    Z3 = Z2; 
    Z2 = Z1; 
    Z1 = nvErro(nIndice);
end; 

// Gera gráfico figure(3);

a3 = figure(3);
a3.background = -2;

plot(Int4p, 'r');

// Gera o grid e as legendas

set(gca(), "background", [-2]); 
xtitle('Integral da forma de onda calculada com 4 pontos', 'Indice das amostras', 'Amplitude'); 
set(gca(), "data_bounds", [0,400,-10,20]); 
set(gca(), "grid", [1,1]); 

.

Os gráficos obtidos são como o mostrado na Figura 5.

.

Os três gráficos obtidos parecem ser iguais, por motivo da forma de onda quadrada ser muito simples. Se olharmos com detalhes os “bicos” das formas de onda ampliadas, poderemos observar as diferenças (Figura 6). Quanto maior o número de pontos usados para o cálculo, menor a diferença se comparada com a onda triangular ideal.

.

Cálculo de Diferencias numéricas utilizando o Scilab

.

Em seguida serão mostradas a utilização das diferenciais calculadas com 2, 3 ou 4 pontos. A onda utilizada para essas operações é a onda triangular resultante da  integral calculada com 4 pontos.

.

// Calculo das diferencias numéricas utilizando-se 2 pontos

// Fecha os gráficos

close(a0);  
close(a1);
close(a2);
close(a3);


Diff2p = zeros(1, size(Int4p,2));
Z1 = 0;
Delta = 1;

// Calcula a Diferencial com 2 pontos

for nIndice = 1:size(nvErro,2)
    Diff2p(nIndice) = ((Int4p(nIndice) - Z1)* 1/Delta);
    Z1 = Int4p(nIndice);
end;

// Apresenta a diferencial no gráfico

a0 = figure(0);
a0.background = -2;

plot(Diff2p, 'm');

// Gera o grid e as legendas

set(gca(), "background", [-2]); 
xtitle('Diferencial da forma de onda calculada com 2 pontos', 'Indice das amostras', 'Amplitude'); 
set(gca(), "data_bounds", [0,400,-1.5,1.5]); 
set(gca(), "grid", [1,1]);

//==========================================================================

// Calculo das diferencias numéricas utilizando-se 3 pontos

Diff3p = zeros(1, size(Int4p,2));
Z1 = 0;
Z2 = 0;
Delta = 1;

// Calcula a Diferencial com 3 pontos

for nIndice = 1:size(nvErro,2)
    Diff3p(nIndice) = ((3*Int4p(nIndice) - 4*Z1 + Z2)* 1/(2*Delta));
    Z2 = Z1;
    Z1 = Int4p(nIndice);
end;

// Apresenta a deferencial no gráfico

a1 = figure(1);
a1.background = -2;

plot(Diff3p, 'b');

// Gera o grid e as legendas

set(gca(), "background", [-2]); 
xtitle('Diferencial da forma de onda calculada com 3 pontos', 'Indice das amostras', 'Amplitude'); 
set(gca(), "data_bounds", [0,400,-2.5,2.5]); 
set(gca(), "grid", [1,1]);

//==========================================================================

// Calculo das diferencias numéricas utilizando-se 4 pontos

Diff4p = zeros(1, size(Int4p,2));
Z1 = 0;
Z2 = 0;
Z3 = 0;
Delta = 1;

// Calcula a Diferencial com 4 pontos

for nIndice = 1:size(nvErro,2)
    Diff4p(nIndice) = ((11*Int4p(nIndice) - 18*Z1 + 9*Z2 - 2*Z3)* 1/(6*Delta));
    Z3 = Z2;
    Z2 = Z1;
    Z1 = Int4p(nIndice);
end;

// Apresenta a deferencial no gráfico

a2 = figure(2);
a2.background = -2;

plot(Diff4p, 'r');

// Gera o grid e as legendas

set(gca(), "background", [-2]); 
xtitle('Diferencial da forma de onda calculada com 4 pontos', 'Indice das amostras', 'Amplitude'); 
set(gca(), "data_bounds", [0,400,-2.5,2.5]); 
set(gca(), "grid", [1,1]);

.

Os gráficos obtidos são semelhantes ao da Figura 7.

.

.

.

Na Figura 8 pode-se observar as três diferenciais sobrepostas. Observe que a diferencial com mais pontos tem um transitório maior do que as outras, o que está coerente. Aqui encerram-se as simulações para demonstrar o funcionamento dos métodos numéricos.

É recomendável que o leitor faça experimentos com o Scilab, crie outras formas de onda, utilize formas de onda reais, pois o Scilab permite que se carregue dados de arquivos, e se familiarize com os recursos do Scilab, se acaso ainda não os conhece. O Scilab é uma ferramenta poderosa, tem recursos para cálculos de problemas de controle, processamento digital de sinais entre muitos outros.

.

Nesse artigo técnico foram desenvolvidos métodos numéricos para o cálculo de integrais e diferenciais, típicas para uso em controle digital por PID. Esses métodos foram simulados no Scilab e foi mostrado o funcionamento correto das diversas aproximações e uma melhor precisão à medida que se utiliza mais pontos para o cálculo dos resultados. As equações desenvolvidas são muito adequadas para a utilização em microprocessadores pela sua simplicidade. Experimente usá-las.

Saiba mais sobre PID Digital

Leia também outros artigos sobre PID já publicados aqui no Embarcados:

• Controle PID em sistemas embarcados, escrito por Carlos Márcio Freitas;
• Controlador PID digital: Uma modelagem prática para microcontroladores – Parte 1,   escrito por Felipe Neves;
• Controlador PID digital, uma modelagem prática para microcontroladores – Parte 2,  escrito  por  Felipe Neves;