Projeto por Emulação

Introdução

Um sistema contendo tanto um sinal contínuo no tempo quanto um sinal discreto é dito um sistema de dados amostrado ¹.

Assumimos que o período de amostragem é fixo. Na prática, sistemas de controle digitais algumas vezes possuem períodos de amostragem variáveis e/ou diferentes períodos de amostragem entre seus ramos de realimentação. Normalmente, a lógica do computador inclui um sinal de clock que fornece um pulso, ou interrupção, a cada $T$ segundos, um conversor A/D envia uma sequencia numérica binária para o computador a cada instante de tempo que a interrupção ocorre. Outra alternativa de implementação, as vezes referenciada como free-running consiste em acessar o conversor A/D depois que cada ciclo de execução do código foi completada. No primeiro caso, o período de amostragem é preciso e fixo. No último caso, o período de amostragem depende do tamanho do código à ser executado, o qual eventuais laços de repetição e desvios condicionais (IFs) devem ser evitados, o que caso contrário, poderia fazer variar a quantidade de código à ser executada, fazendo variar o período de amostragem.

Também num sistema de controle digital deve estar presente um sampler (amostrador) e um conversor A/D para a entrada de comando (ou referência), $r(t)$, o que gera o sinal discretizado $r(kT)$ que juntamente com o saída amostrada do sistema, $y(kT)$ gera o sinal de erro discretizado: $e(kT)=r(kT)-y(kT)$:

Um sistema contínuo de compensação pode ser aproximado por uma equação de diferenças que seria uma versão discreta de uma equação diferencial e pode ser adotado para emular o comportamento dinâmico de um compensador contínuo, $C(s)$, se o período de amostragem for curto o suficiente. O resultado desta equação de diferenças seria um sinal discreto $u(kT)$ a cada instante de amostragem. Este sinal seria convertido para o mundo contínuo, $u(t)$ através de um conversor D/A e um sustentador (Holder). O conversor D/A converte o número binário numa tensão analógica e o sustentador manteria a última tensão analógica gerada durante os instantes de amostragem.

Existem duas técnicas básicas para determinar a equação de diferenças que define o controlador digital. Uma delas, chamada emulação, consistem em projetar o compensador contínuo no tempo, $C(s)$ usando métodos clássicos de controle automático e então realizar uma aproximação de $C(s)$ para uma representação discreta, $C(z)$ usando um dos seguintes métodos:

• Método da Correspondência Pólo-Zero (MPZ = Matched Pole-Zero).
• Método de Euler;
• Método de Tustin;

Método da Correspondência Pólo-Zero (ou MPZ)

Este método, o da correspondência pólo-zero leva em conta a definição da transformada $\mathcal{Z}$, extrapolando então a relação entre os planos $s$ e $z$, obedecendo a relação:

$z=e^{sT} \quad \quad$ (Pela definição)

como: $s=\alpha+j\omega$, então teremos:

$z=e^{T(\alpha+j\omega)}=e^{\alpha T} \, e^{j \omega T}$

$z=e^{\alpha T} ( \cos{\omega T} + j \sin{\omega T} )$

$z=e^{\alpha T} \, \angle \, {\omega T}$

uma vez que: $( \cos{\omega T} + j \sin{\omega T} ) = 1 \, \angle \, {\omega T}$.

Desta forma, cada região do plano-s pode ser mapeada na correspondente região do plano-z.

Se tomamos a transformada $\mathcal{Z}$ da função amostrada $x(k)$, então os pólos de $X(z)$ estão relacionados com os pólos de $X(s)$ obedecendo a relação $z=e^{sT}$. O problema é que temos que considerar também que $s=(1/T)\ln{z}$ é a inversa de $z=e^{sT}$. Sendo assim, se temos $P(z)$ e substituímos $z=e^{Ts}$, obtemos $P(e^{Ts})$ como resultado da conversão para o plano$-s$. De modo simular, se temos $P(s)$ e substituímos $s=(1/T)\ln{z}$, obtemos $P[(1/T)\ln{z}]$ como resultado da conversão para o plano$-z$. Infelizmente, ambas transformações podem ser levadas à caso, respeitando detalhes da complicada transformada mas conduzem a funções trancedentais (Cap. 13 de ²).

Relações entre o plano-s e o plano-z

A figura à seguir mostra as linhas de mapeamento resultantes para valores constantes de $\zeta$ e para valores constantes de $w_n$ baseado na relação: $z=e^{sT}$:

E a próxima figura mostra diferentes respostas no tempo resultantes para diferentes posições para um pólo em malha fechada.

Método de Euler

Uma maneira particularmente simples de fazer com que um computador digital se aproxime da solução em tempo-real de uma equação diferencial é usar o método de Euler (também conhecido como regra retangular direta -- forward rectangular rule). Ela é resultado da seguinte aproximação:

$\dfrac{dx}{dt}=\displaystyle\lim_{\Delta t \to 0}{\dfrac{\Delta x}{\Delta t}}$

que aplicada neste caso em particular resulta em:

$\dot{x} \approx \dfrac{x(k)-x(k-1)}{T}$

onde:
$T$: corresponde ao período de amostragem (em segundos), $T = \Delta t = t(k)-t(k-1)$;
$x(k)$: valor de $x$ no instante $t_k$, ou $x(kT)$;
$x(k-1)$: valor de $x$ no instante $t_{k-1}$, ou $x[(k-1)T]$.

Esta aproximação pode ser usada no lugar de várias derivadas que apareçam na equação diferencial do controlador e resulta num conjunto de equações que pode ser calculada por um computador digital. Estas equações de diferenças são resolvidas repetidamente a cada $T$ período de tempo. Para sistemas com margem de banda da ordem de uns poucos Hertz, taxas de amostragem resultam normalmente na ordem de 1 KHz, e assim os períodos de amostragem ficarão em torno de 1 ms e os erros envolvidos na aproximação resultarão bastante pequenos ¹.

Exemplo_1: Caso de um Controlador de 1a-ordem

Usando o método de Euler, encontre a equação de diferenças que necessita ser programada num computador se $C(s)$ é dado como:

$C(s)=\dfrac{U(s)}{E(s)}=K_o \, \dfrac{s+a}{s+b}$

Solução: Primeiro determinamos a equação diferencial que corresponde à $C(s)$:

$(s+b)U(s)=K_o (s+a)E(s)$

por inspeção percebemos que a correspondente equação diferencial é:

$\dot{u}+bu=K_o (\dot{e}+ae) \quad \quad \text{(1)}$

Usando o método de Euler para aproximar a eq. (1) obtermos a seguinte equação de diferenças:

$\dfrac{u(k)-u(k-1)}{T}+bu(k)=K_o \left[ \dfrac{e(k)-e(k-1)}{T} + ae(k) \right] \quad \quad \text{(2)}$

Rearrumando a eq. (2) obtemos:

$u(k)=u(k-1)+T\left\{ (-b) \cdot u(k) + K_o \left[ \dfrac{e(k)-e(k-1)}{T} + a \cdot e(k)\right] \right\} \quad \quad \text{(3)}$

A equação (3) determina o novo valor de controle, $u(k)$, uma vez conhecida a amostra passada do sinal de controle, $u(k-1)$ e os novos e valores passados do sinal de erro, $e(k)$ e $e(k-1)$.

Exemplo_2: Caso de um Controlador por Avanço de Fase

Encontre uma implementação digital para compensador de avanço de fase definido para a seguinte planta (Exemplo 5.13 de ¹):

$G(s)=\dfrac{1}{s(s+1)}$

usando taxas de amostragem de: a) 20 e b) 40 Hz, simular as equações de controle para uma entrada degrau e comparar o resultado obtido entre o controlador analógico e o digital.

O sistema deve responder em malha-fechada com fator de amortecimento $\zeta>0,5$ e com frequência natural $w_n>7$ (rad/s).

Na ocasião deste exercício foram feitas diferentes tentativas de pares pólo-zero para o compensador de avanço de fase como mostra a próxima figura:

(a) $C(s)=(s+2)/(s+20)$	(b) $C(s)=(s+2)/(s+10)$
(c) $C(s)=s+2$

Na ocasião foi escolhido o controlador por atraso de fase:

$C(s)=\dfrac{70(s+2)}{(s+10}$

A função transferência em malha aberta resulta:

$FTMA(s)=C(s)G(s)=\dfrac{70(s+2)}{s(s+10)(s+1)}$

onde $K=70$. Este ganho resulta na resposta ao degrau mostrada na seguinte figura:

Solução: a função transferência do compensador de avanço de fase para esta planta fica (baseado no exemplo 5.13 de ¹):

$C(s)=\dfrac{70(s+2)}{s+10}$

Desta forma, os parâmetros para a eq. (3) são $a=2$, $b=10$ e $K=70$.

a) Para a taxa de amostragem de 20 Hz, $T=0,05$ (segundos) e a eq. (3) pode ser simplificada para:

$u(k)=0,5u(k-1)+70[e(k)-0,9e(k-1)]$

Colocando no formato de $C(z)$ temos:

$U(z)=0,5z^{-1}U(z)+70[E(z)-0,9z^{-1}E(z)]$

trabalhando a expressão anterior obtemos:

$U(z)[1-0,5z^{-1}]=70E(z)[1-0,9z^{-1}]$

e então:

$C(z)=\dfrac{U(z)}{E(z)}=\dfrac{70(1-0,9z^{-1})}{(1-0,5z^{-1})}=\dfrac{70(z-0,9)}{(z-0,5)}$

b) Para a taxa de amostragem de 40 Hz, $T=0,025$ (segundos), a eq. (3) simplificada resulta em:

$u(k)=0,75u(k-1)+70[e(k)-0,95e(k-1)]$

Colocando no formato de $C(z)$ obtemos:

$C(z)=\dfrac{U(z)}{E(z)}=\dfrac{70(1-0,95z^{-1})}{[1-0,75z^{-1}]}=\dfrac{70(z-0,95)}{(z-0,75)}$

O código no MATLAB para calcular a resposta contínua ao degrau fica:

xxxxxxxxxx
>> num=70*[1 2]
>> den=conv([1 1 0],[1 10])
>> [numcl,dencl]=feedback(num,den,1,1)
>> step (numcl,dencl)

Que resulta no gráfico mostrado na seguinte figura anterior.

A próxima figura mostra o diagrama em blocos do sistema simulado usando tanto o compensador contínuo no tempo (a) quanto o simulador discreto (b).

As respostas ao degrau para (a) T=0.05 e (b) T=0.025 aparecem na próxima figura:

(a) $f_s=20$ [Hz]	(b) $f_s=40$ [Hz]

Conclusão: Este método (de Euler) é mais adequado quanto menor for o período de amostragem, $T$. Para grandes valores de $T$, este método pode resultar em instabilidade, mesmo quando o sistema sendo discretizado é estável no tempo contínuo.

Método de Tustin

Outro método de discretização é tratar o problema como uma integração numérica (aproximação).

Suponha que:

$\dfrac{U(s)}{E(s)}=C(s)=\dfrac{1}{s}$

$o qual é uma integração. Assim sendo:

$u(kT)=\displaystyle\int_{0}^{kT-T}{e(t)dt}+\displaystyle\int_{kT-T}^{kT}{e(t)dt}$

pode ser reescrito como:

$u(kT)=u(kT-T)+$ área entre $e(t)$ e período de amostragem atual $kT \quad \quad \text{(4)}$

O método de Tustin realiza integração numérica trapezoidal a cada passo do período de amostragem, tentando aproximar $e(t)$ através de uma linha direta entre dois instantes de amostragem:

Ver PID Digital: Integração Trapezoidal, que resulta na mesma equação mostrada abaixo.

Reescrevendo $u(kT)$ como $u(k)$ e $u(kT-T)$ como $u(k-1)$ para simplificar a escrita, tra nsformamos a equação (4) em:

$u(k)=u(k-1)+\dfrac{T}{2}\left[ e(k-1)+e(k) \right]$

tomando a transformada$-\mathcal{Z}$ da expressão anterior, obtemos:

$\dfrac{U(z)}{E(z)}=\dfrac{T}{2}\left( \dfrac{1+z^{-1}}{1-z^{-1}} \right)=\dfrac{1}{ \dfrac{2}{T}\left( \dfrac{1-z^{-1}}{1+z^{-1}} \right) }$

Se: $C(s)=\dfrac{a}{(s+a)}$, aplicando esta aproximação para a integração resulta em:

$C(z)=\dfrac{a}{ \dfrac{2}{T} \left( \dfrac{1-z^{-1}}{1+z^{-1}} \right) + a }$

Note que se substituírmos:

$s=\dfrac{2}{T}\left( \dfrac{1-z^{-1}}{1+z^{-1}} \right) \quad \quad \text{(5)}$

para cada ocorrência de $s$ em qualquer $C(s)$, obtemos $C(z)$ baseado na equação da integração trapezoidal.

Este é o chamado método de transformação bilinear de Tustin ¹. A idéia é a de que nos instantes amostrados, a transformação bilinear seja capaz de preservar o mesmo tipo (valor) da resposta de um controlador contínuo no tempo. Em resumo este tipo de transformação é usado para transformar um compensador contínuo no tempo, $C(s)$ na sua versão digital, $C(z)$. Explicitamente a transformação bilinear de Tustin é dada por:

$s=\dfrac{2}{T} \, \dfrac{(z-1)}{(z+1)} \quad \quad \text{(6)}$

e sua inversa é dada por:

$\begin{array}{rcl}z &=& \dfrac{-\left( s+ \dfrac{2}{T} \right)}{\left( s - \dfrac{2}{T} \right)}\\ &=& \dfrac{1+\dfrac{T}{2}s}{1-\dfrac{T}{2}s} \end{array} \quad \quad \text{(7)}$

Quanto menor for o intervalo de amostragem, $T$ (ou maior a frequência de amostragem), mais a versão discretizada se aproxima do comportamento no domínio-frequência do compensador no domínio contínuo de tempo. A eq. (6) é a que entre os métodos de integração numérica, permite levar a um resultado mais próximo da resposta (em tempo e frequência) que realmente seria obtido usando o compensador no domínio contínuo de tempo.

Åström e Wittenmark³, propuseram uma referência para selecionar o período de amostragem à ser adotado. Eles concluíram que o valor de $T$ em segundos deve estar na faixa de $0,15/\omega_{\Phi M}$ até $0,5/\omega_{\Phi M}$, onde $\omega_{\Phi M}$ corresponde a frequência (em rad/s) na qual o diagrama de magnitude da função transferência analógica do compensador mais planta passa por 0 dB (ganho unitário).

Exemplo_1: Equação de diferenças usando o Método de Tustin

Voltando ao exemplo da seção anterior e usando taxa de amostragem de $f_s=15$ [Hz], podemos usar o MATLAB para computar a aproximação de Tustin.

Solução: A forma contínua do compensador estudado na seção anterior era:

$C(s)=70 \, \dfrac{(s+2)}{(s+10)}$

Usando declarações na linha de comandos do MATLAB, podemos fazer:

xxxxxxxxxx
>> numc=70*[1 2];
>> denc=[1 10];
>> T=1/15;
>> C=tf(numc,denc);
>> zpk(C)
Zero/pole/gain:
70 (s+2)
--------
 (s+10)
>> Cd=c2d(C,T,'tustin') 
Transfer function:
56 z - 49
---------
 z - 0.5
Sampling time: 0.066667
>> zpk(Cd)
Zero/pole/gain:
56 (z-0.875)
------------
  (z-0.5)
Sampling time: 0.066667
>> 

que produz:

$C(z)=\dfrac{56-49z^{-1}}{1-0,5z^{-1}}=\dfrac{56 (z-0.875)}{(z-0.5)}$

que pode ser transcrito para uma equação de diferenças resultando em:

$u(k)=0,223u(k-1)+70e(k)-59,12e(k-1)$

que é similar a equação obtida usando o método de Euler.

Problema sugerido Comparar a resposta deste controlador usando o método de Tustin com os resultados obtidos quando o mesmo compensador contínuo no tempo foi discretizado usando o método de Euler (seção anterior). Para os mesmos 2 valores de período de amostragem: a) $f_s=20$ (Hz) e b) $f_s=40$ (Hz).

Simulando os controladores no Simulink, arquivo: Lead_Tustin_2a_ordem.slx:

Resultados obtidos:

Exemplo_2: Projeto de Controlador Digital usando o Método de Tustin

Seja um sistema (planta) descrito pela equação (Exemplo 13.12, Cap 13 de ²):

$G(s)=\dfrac{1}{s(s+6)(s+10)} \quad \quad \text{(8)}$

onde o projeto de um compensador por avanço resultou num sistema compensado que opera com sobre sinal máximo de 20% e um tempo de acomodação de 1,1 segundos. Defina a equação de diferenças que corresponde ao compensador contínuo transformado para digital usando o método de Tustin. A equação do compensador é dada por:

$C(s)=\dfrac{1977(s+6)}{(s+29,1)} \quad \quad \text{(9)}$

Solução: Com base na eq. (8) e (8) se percebe que a frequência que corresponde à 0 dB, $\omega_{\Phi M}$, para $C(s)G(s)$, é de 5,8 rad/s. Com base na linha guia proposta por ³, o menor valor recomendado para $T$ deveria estar na faixa de $0,15/\omega_{\Phi M}=0,026$ até $0,5/\omega_{\Phi M}=0,086$ segundos. Iremos adotar o valor de $T=0,01$ segundos.

Substituindo a eq. (6) na eq. (9) com $T=0,01$ segundos, leva à:

$C(z)=\dfrac{1778z-1674}{z-0,746}=\dfrac{1777.7 (z-0.9417)}{(z-0.746)} \quad \quad \text{(10)}$

e a função transferência da planta usando sustentador de ordem zero (ZOH), com $T=0,01$ segundos resulta em:

$BoG(z)=\dfrac{(1,602 \times 10^{-7}z^2)+(6,156\times 10^{-7}z)(1,478 \times 10^{-7})}{z^3-2,847z^2+2,699z-0,8521} \quad \quad \text{(11)}$

Obs.: As equações anteriores foram obtidas usando o seguinte código no MATLAB:

xxxxxxxxxx
>> num=1;
>> den=poly([0 -6 -10]);
>> G=tf(num,den);
>> zpk(G)
        1
  --------------
  s (s+10) (s+6) 
Continuous-time zero/pole/gain model.
>> num_c=1977*[1 6];
>> den_c=[1 29.1];
>> C=tf(num_c,den_c);
>> zpk(C)
  1977 (s+6)
  ----------
   (s+29.1) 
Continuous-time zero/pole/gain model.
>> % convertendo para controle digital
>> T=0.01;
>> help c2d
 c2d  Converts continuous-time dynamic system to discrete time.
 
    SYSD = c2d(SYSC,TS,METHOD) computes a discrete-time model SYSD with 
    sampling time TS that approximates the continuous-time model SYSC.
    The string METHOD selects the discretization method among the following:
       'zoh'       Zero-order hold on the inputs
       'foh'       Linear interpolation of inputs
       'impulse'   Impulse-invariant discretization
       'tustin'    Bilinear (Tustin) approximation.
       'matched'   Matched pole-zero method (for SISO systems only).
    The default is 'zoh' when METHOD is omitted. The sampling time TS should 
    be specified in the time units of SYSC (see "TimeUnit" property).
​
>> Cd=c2d(C,T,'tustin');
>> zpk(Cd)
  1777.7 (z-0.9417)
  -----------------
      (z-0.746)
Sample time: 0.01 seconds
Discrete-time zero/pole/gain model.
>> % convertendo a planta para o mundo digital usando ZOH:
>> Gd=c2d(G,T)
Gd =
  1.602e-07 z^2 + 6.156e-07 z + 1.478e-07
  ---------------------------------------
    z^3 - 2.847 z^2 + 2.699 z - 0.8521
Sample time: 0.01 seconds
Discrete-time transfer function.
>> zpk(Gd)
  1.6016e-07 (z+3.586) (z+0.2574)
  -------------------------------
    (z-1) (z-0.9418) (z-0.9048)
Sample time: 0.01 seconds
Discrete-time zero/pole/gain model.
>> ftma=series(C,G);        % FTMA(s)
>> ftmad=serires(Cd,Gd);    % FTMA(z) com T=0,01 (s)
>> ftmf=feedback(ftma,1);   % FTMF(s)
>> ftmfd=feedback(ftmad,1); % FTMF(z)
>> [y,t]=step(ftmf);        % resposta degrau, sistema contínuo
>> plot(t,y)
>> [yd1,td1]=step(ftmfd);
>> hold on;
>> stairs(td1,yd1,'m-')
>> legend('FTMF(s)','FTMF(z) com T=0,01')
>> title('Resposta ao Degrau');
>> xlabel('Tempo (s)');
>> ylabel('y^*(t)')
>> grid

A resposta para entrada degrau para o sistema compensando é mostrada na figura à seguir, que compara as respostas do compensador contínuo no tempo versus compensadores digitais usando $T=0,01$ segundos.

(a) Resposta ao degrau.	(b) "Zoom'' sobre o gráfico (a).

Outras Transformações Bilineares

Note que o método de Tustin trata-se de um caso específico da aplicação de um método de transformação bilinear do tipo:

$z=\dfrac{as+b}{cs+d}$

e sua inversa:

$s=\dfrac{-dz+b}{cz-a}$

A idéia é sempre buscar uma simples transformação capaz de atuar em ambos os lados (plano-s e plano-z) através de algum método de substituição direta (Cap 13. de ²).

Por exemplo, a transformação bilinear:

$z=\dfrac{s+1}{s-1} \qquad \mbox{e sua inversa: } \qquad s=\dfrac{z+1}{z-1}$

é proposta por Kuo⁴ para aproximar raízes presentes no interior do círculo unitário $|z|=1$ para o semi-plano esquerdo do plano-s. No caso, esta transformação recebeu o nome de transformação-R.

Outro tipo de transformação bilinear é a transformada-w onde:

$\omega=\dfrac{1}{\ln{(z)}}$

tenta transformar o círculo unitário do plano-z para o eixo imaginário de outro plano complexo $\omega$ ⁵. Neste caso:

$\omega=\dfrac{2}{T}\dfrac{(z-1)}{(z+1)} \quad \quad \text{(12)}$

e sua inversa:

$\begin{array}{rcl}z&=&\dfrac{ \dfrac{2}{T}+\omega }{ \dfrac{2}{T}-\omega }\\ &=&\dfrac{1+\dfrac{\omega \, T}{2}}{1-\dfrac{\omega \, T}{2}} \end{array} \quad \quad \text{(13)}$

que é uma variação da proposta realizada pelo Ogata (Cap. 4: Design of Discrete-Time Control Systems by Conventional Methods, de ⁶):

$z=\dfrac{w+1}{w-1} \qquad \text{e sua inversa: } \qquad w=\dfrac{z+1}{z-1}$

Em comum as últimas 2 transformadas permitem realizar todo um projeto de controlador discreto usando o mapeamento em frequencia, tal qual se fazia usando Diagramas de Bode quando estávamos no plano-s (ou no mundo contínuo). Neste último caso, estas transformadas permitem converter $BoG(z)$ em $BoG(w)$ e lidar com o projeto do controlador tal qual se fazia antes usando Diagrama de Bode -- o que muda agora é que ao invés do Diagrama de Bode ser feito no plano-$s$ é realizado sobre o plano-$w$.

Para "imprimir" figuras de modelos do Simulink

Você pode usar a linha de comandos para transformar numa figura, um modelo do Simulink sem recorrer a um print screen com resultados deploráveis (baixa resolução, 75 dpi). Para tanto:

xxxxxxxxxx
>> print('-sLead_Tustin_2a_ordem', '-r150', '-dpng','Lead_Tustin_2a_ordem.png')  % gera arquivo .png
>> print('-sLead_Tustin_2a_ordem', '-r150', '-djpeg','Lead_Tustin_2a_ordem.jpg') % gera arquivo .jpg
>> print('-sLead_Tustin_2a_ordem', '-r150', '-dpdf','Lead_Tustin_2a_ordem.pdf')  % gera arquivo .pdf

Detalhes:

• print(-s<Model_Name>,...) indica que se deseja "imprimir" não uma janela gráfica, mas sim um modelo do Simulink;
• a opção '-r150', indica a resolução desejada, neste caso, 150 dpi;
• o último parâmetro, se refere ao nome do arquivo.

Neste caso, foram gerados:

xxxxxxxxxx
-rw-r--r--  1 fernandopassold  staff  79721  2 Jun 16:42 Lead_Tustin_2a_ordem.jpg
-rw-r--r--  1 fernandopassold  staff  28346  2 Jun 16:43 Lead_Tustin_2a_ordem.pdf
-rw-r--r--  1 fernandopassold  staff  66775  2 Jun 16:42 Lead_Tustin_2a_ordem.png

Fernando Passold (primeira publicação: 02/06/2014; atualizado em 02/06/2025)