Projeto direto de controladores digitais

PI, PI Tipo 2 e PID a partir do modelo discreto G(z)

Prof. Dr. Raphael Teixeira

UFPA — Campus Tucuruí — Faculdade de Engenharia Elétrica

Roteiro

Emulação × projeto direto
Plantas de teste e especificações no plano \(z\)
Método 1 — Alocação polinomial (equação diofantina): PI e PID
Método 2 — Lugar das raízes em \(z\): PI
Método 3 — Resposta em frequência no plano \(w\): PI Tipo 2
Robustez, estrutura 2 GDL e quadro comparativo

Motivação

Emulação × projeto direto

Emulação

Projeta \(C(s)\) sobre \(G(s)\) e discretiza (Tustin, ZOH, …)
Exige um modelo contínuo
Aproximação válida só para \(\omega_c \ll \omega_s\)
Atraso do ZOH e computacional precisam ser adicionados “à mão”

Projeto direto

Projeta \(C(z)\) sobre \(G(z)\)
O ARX já é um \(G(z)\): nada de d2c
ZOH, atraso de transporte \(n_k\) e atraso computacional já estão no modelo
Resultado exato no instante de amostragem

Importante

Converter o ARX para \(s\) com d2c pode ser mal condicionado (polos reais negativos não têm correspondente contínuo) e reintroduz aproximações.

Do ARX ao \(G(z)\)

\[ A(q^{-1})\,y(k) = B(q^{-1})\,u(k-n_k) + e(k) \quad\Longrightarrow\quad G(z) = \frac{B(z)}{A(z)} \]

% dados: objeto iddata (entrada PRBS, saída medida), Ts = T
m  = arx(dados, [na nb nk]);
Gz = tf(m);                     % apenas o canal medido (sem ruído)
[B, A] = tfdata(Gz, 'v');       % num/den com mesmo comprimento,
                                % potências decrescentes de z

Nesta aula, as plantas são obtidas por ZOH exato. Sem ruído, o ARX com \([n_a\ n_b\ n_k]\) corretos recupera exatamente esses modelos.

Plantas e especificações

Plantas de teste (\(T = 10\) ms)

Planta 1 — 1ª ordem \[ G_1(s) = \frac{2}{0{,}1s+1} \] \[ G_1(z) = \frac{0{,}1903}{z-0{,}9048} \]

Planta 2 — 2ª ordem (\(\omega_n = 10\), \(\zeta = 0{,}5\)) \[ G_2(s) = \frac{100}{s^2+10s+100} \] \[ G_2(z) = \frac{0{,}004833\,z + 0{,}004675}{z^2 - 1{,}8953\,z + 0{,}9048} \]

T  = 0.01;  s = tf('s');
G1 = 2/(0.1*s + 1);            G1z = c2d(G1, T, 'zoh');
G2 = 100/(s^2 + 10*s + 100);   G2z = c2d(G2, T, 'zoh');
[B1, A1] = tfdata(G1z, 'v');
[B2, A2] = tfdata(G2z, 'v');

Da especificação ao plano \(z\)

Polos contínuos desejados \(\;\longrightarrow\;\) mapeamento exato \(z = e^{sT}\):

\[ s_d = -\zeta\omega_n \pm j\omega_n\sqrt{1-\zeta^2} \qquad\Rightarrow\qquad z_d = e^{s_d T} \]

\[ P_d(z) = (z - z_d)(z - \bar z_d)\,\prod_i (z - z_{a,i}) \]

Os polos auxiliares \(z_{a,i}\) devem ser mais rápidos que o par dominante; aqui, \(z_a = e^{-4\zeta\omega_n T}\).

polo = @(zeta, wn) exp(T*(-zeta*wn + 1j*wn*sqrt(1 - zeta^2)));

Graus de liberdade: quem resolve o quê

Controlador com integrador: \(C(z) = \dfrac{S(z)}{(z-1)\,R'(z)}\)

Planta	Controlador	ordem de \(P_d\)	incógnitas	Consequência
1ª ordem	PI	2	2 (\(s_0, s_1\))	solução exata
1ª ordem	PID (\(r_1 = 0\))	3	3	exata, mas \(K_d \le 0\)
2ª ordem	PI	3	2	subdeterminado → LGR / frequência
2ª ordem	PID (\(r_1\) livre)	4	4	solução exata
qualquer	Tipo 2	—	\(K_c, \omega_z, \omega_p\)	frequência (fator \(k\))

Regra geral: a solução é única quando \(\deg S = \deg A\) e \(\deg R' = \deg B\), já contando o integrador.

Método 1 — Alocação polinomial

A equação diofantina

\[ \underbrace{A(z)\,(z-1)\,R'(z)}_{A(z)R(z)} + B(z)\,S(z) = P_d(z) \]

\((z-1)\): integrador imposto, garante erro nulo ao degrau
\(R'(z)\) mônico: polos livres do controlador (ex.: filtro da derivada)
Igualando coeficientes, obtém-se um sistema linear (matriz de Sylvester)

A lei de controle na forma RST, com \(T(z)\) a ser discutido adiante:

\[ R(q)\,u(k) = T(q)\,r(k) - S(q)\,y(k) \]

Função `rst_poly.m`

function [R, S] = rst_poly(A, B, Rfix, nRp, nS, P)
% Resolve  A*Rfix*R' + B*S = P
%   Rfix : parte fixa de R (ex.: [1 -1] = integrador)
%   nRp  : grau de R' (mônico, coeficientes livres)
%   nS   : grau de S
%   Polinômios em potências decrescentes de z
Ab  = conv(A, Rfix);
n   = numel(P) - 1;
pad = @(p) [zeros(1, n + 1 - numel(p)) p];
M = [];
for i = 1:nRp
    M = [M, pad(conv(Ab, [1 zeros(1, nRp - i)])).'];
end
for j = 0:nS
    M = [M, pad(conv(B,  [1 zeros(1, nS - j)])).'];
end
rhs = (P - pad(conv(Ab, [1 zeros(1, nRp)]))).';
x   = M(2:end, :) \ rhs(2:end);      % 1ª linha: 1 = 1 (P mônico)
R   = conv(Rfix, [1 x(1:nRp).']);
S   = x(nRp + 1:end).';
end

Caso 1 — PI na planta 1 (solução analítica)

\(G_1(z) = \dfrac{b}{z-a}\), \(\quad C(z) = \dfrac{s_0 z + s_1}{z-1}\)

\[ (z-1)(z-a) + b(s_0 z + s_1) = z^2 + p_1 z + p_2 \;\Rightarrow\; s_0 = \frac{p_1 + 1 + a}{b}, \quad s_1 = \frac{p_2 - a}{b} \]

Na forma paralela (Euler para trás), \(C = K_p + K_i T\dfrac{z}{z-1}\), o que dá \(K_p = -s_1\) e \(K_i = (s_0+s_1)/T\).

Especificação: \(\zeta = 0{,}8\) e \(t_s = 0{,}15\) s, logo \(\omega_n = 33{,}3\) rad/s e \(z_d = 0{,}7507 \pm j\,0{,}1522\).

zd = polo(0.8, 4/(0.8*0.15));
Pd = real(poly([zd conj(zd)]));             % z^2 - 1.5013 z + 0.5867
[R, S] = rst_poly(A1, B1, [1 -1], 0, 1, Pd);
Kp = -S(2);   Ki = sum(S)/T;

Resultado: \(s_0 = 2{,}1201\) e \(s_1 = -1{,}6718\), ou seja, \(K_p = 1{,}672\) e \(K_i = 44{,}83\).

Caso 2 — PID na planta 1: o que a álgebra diz

Com \(R = (z-1)\,z\) (derivada sem filtro) e \(P_d\) de 3ª ordem:

\[ s_0 = \frac{p_1 + 1 + a}{b}, \qquad s_1 = \frac{p_2 - a}{b}, \qquad s_2 = \frac{p_3}{b} \]

Os ganhos seguem de \(s_0 = K_p + K_iT + K_d/T\), \(\;s_1 = -K_p - 2K_d/T\) e \(\;s_2 = K_d/T\).

\(p_3 = -|z_d|^2 z_a\). Para um polo auxiliar real e positivo, \(p_3 < 0\), então \(K_d < 0\).
Com \(z_a = 0{,}344\), obtém-se \(K_d = -0{,}0106\).
Com \(z_a = 0\), obtém-se \(K_d = 0\), e o controlador volta a ser o PI.

Nota

Conclusão: para uma planta de 1ª ordem, a ação derivativa é redundante. A contagem de graus de liberdade mostra isso antes de qualquer simulação.

Caso 3 — PID na planta 2

\[ C(z) = \frac{s_0 z^2 + s_1 z + s_2}{(z-1)(z-p)}, \qquad p:\ \text{polo do filtro da derivada (incógnita)} \]

Especificação: \(\zeta = 0{,}7\) e \(\omega_n = 15\) rad/s, o que dá \(z_d = 0{,}8952 \pm j\,0{,}0963\). Os polos auxiliares ficam em \(z_a = 0{,}657\) (duplo).

zd = polo(0.7, 15);   za = exp(-4*0.7*15*T);
P4 = real(poly([zd conj(zd) za za]));
[R, S] = rst_poly(A2, B2, [1 -1], 1, 2, P4);
p = R(3);             % R = z^2 - (1+p) z + p

Resultado: \(S = [\,19{,}304\;\; -35{,}378\;\; 16{,}325\,]\) e \(R = (z-1)(z-0{,}3024)\).

Aqui \(p > 0\), o que corresponde a um filtro fisicamente coerente. Se \(p < 0\) (polo no semieixo real negativo), \(u(k)\) oscila a cada amostra; nesse caso, reposicione os polos auxiliares.

Conversão para ganhos paralelos

Com integral e derivada filtrada por Euler para trás:

\[ C(z) = K_p + K_i T\frac{z}{z-1} + \frac{K_d}{T_f + T}\,\frac{z-1}{z-p}, \qquad p = \frac{T_f}{T_f + T} \]

Igualando os numeradores sobre \((z-1)(z-p)\), com \(K_d' = K_d/(T_f+T)\):

\[ \begin{bmatrix} 1 & T & 1 \\ -(1+p) & -pT & -2 \\ p & 0 & 1 \end{bmatrix} \begin{bmatrix} K_p \\ K_i \\ K_d' \end{bmatrix} = \begin{bmatrix} s_0 \\ s_1 \\ s_2 \end{bmatrix} \]

g  = [1 T 1; -(1+p) -p*T -2; p 0 1] \ S.';
Kp = g(1);  Ki = g(2);  Kd = g(3)*T/(1-p);  Tf = p*T/(1-p);
Cp = pid(Kp, Ki, Kd, Tf, T, 'IFormula','BackwardEuler', ...
                             'DFormula','BackwardEuler');
tf(Cp) - tf(S, R, T)          % deve ser ~0 (verificação)

Resultado: \(K_p = 3{,}756\), \(K_i = 35{,}92\), \(K_d = 0{,}2177\) e \(T_f = 4{,}33\) ms.

Zeros de \(S\) e a estrutura 2 GDL

Com \(T(z) = S(z)\) (1 GDL, erro na entrada do controlador):

\[ \frac{Y}{R} = \frac{B\,S}{P_d} \quad\longrightarrow\quad \text{zeros de } S \text{ causam sobressinal} \]

Com \(T(z) = t_0 = S(1)\) (2 GDL):

\[ \frac{Y}{R} = \frac{t_0\,B}{P_d}, \qquad \left.\frac{Y}{R}\right|_{z=1} = 1 \quad (\text{pois } R(1) = 0) \]

No PI, isso equivale à estrutura I-P: a referência entra só pelo termo integral.
A dinâmica de rejeição de perturbação não muda, porque ela depende apenas de \(P_d\).

Implementação no DSP (forma RST)

PID da planta 2, com \(R = z^2 - (1+p)z + p\):

\[ \begin{aligned} u(k) ={}& (1+p)\,u(k-1) - p\,u(k-2) \\ & + t_0\,r(k) - s_0\,y(k) - s_1\,y(k-1) - s_2\,y(k-2) \end{aligned} \]

com \(t_0 = s_0 + s_1 + s_2 = 0{,}2506\).

A forma incremental facilita o anti-windup: satura-se \(u(k)\) e armazena-se o valor saturado em \(u(k-1)\).
Os coeficientes saem diretamente da diofantina, sem passar pelos ganhos paralelos.

Método 2 — Lugar das raízes em \(z\)

Caso 4 — PI na planta 2

A diofantina é subdeterminada (3 polos, 2 incógnitas). O procedimento é fixar o zero do PI e escolher o ganho pelo LGR.

\[ C(z) = K\,\frac{z - z_c}{z - 1}, \qquad z_c = e^{-\omega_i T}, \quad \omega_i = 8\ \text{rad/s} \Rightarrow z_c = 0{,}9231 \]

zc = exp(-8*T);
L0 = tf([1 -zc], [1 -1], T) * G2z;
figure; rlocus(L0); zgrid; axis([0.85 1.02 -0.2 0.2])

% ganho: menor raio dos polos em MF, com zeta >= 0.4 em todos
best = [Inf NaN];
for K = linspace(0.01, 2, 4000)
    pc = pole(feedback(K*L0, 1));
    sc = log(pc)/T;  zet = -real(sc)./abs(sc);
    if max(abs(pc)) < 1 && min(zet) >= 0.4 && max(abs(pc)) < best(1)
        best = [max(abs(pc)) K];
    end
end
K = best(2);   S = K*[1 -zc];   R = [1 -1];

Resultado: \(K \approx 0{,}224\).

Leitura do resultado

\(t_s \approx 2{,}4\) s em malha fechada, contra \(\approx 0{,}8\) s em malha aberta: o PI deixa a planta mais lenta.
O polo do integrador fica preso perto de \(z = 1\), porque o par complexo da planta limita o ganho.
A robustez é excelente (\(M_s = 1{,}20\)), mas o desempenho é pobre.

Dica

O LGR mostra claramente o limite estrutural do PI: sem um segundo zero, não há como puxar o par complexo para a esquerda. Esse é o argumento natural para o PID, ou para o Tipo 3.

Método 3 — Resposta em frequência no plano \(w\)

A transformação \(w\)

\[ w = \frac{2}{T}\,\frac{z-1}{z+1} \qquad\Longleftrightarrow\qquad z = \frac{1 + wT/2}{1 - wT/2} \]

Sobre o círculo unitário, \(w = j\nu\) com \(\nu = \dfrac{2}{T}\tan\!\left(\dfrac{\omega T}{2}\right)\).
\(G_w(j\nu) = G(e^{j\omega T})\) exatamente: o Bode no plano \(w\) é o Bode do modelo discreto.
O ZOH e o atraso aparecem como um zero de fase não mínima em \(w = 2/T\).
Em MATLAB: Gw = d2c(Gz,'tustin'). Projeta-se \(C_w\) com as regras clássicas e volta-se com c2d(Cw,T,'tustin').

PI Tipo 2 e o fator \(k\) (Venable)

\[ C_w(w) = K_c\,\frac{1 + w/\omega_z}{w\,(1 + w/\omega_p)} \]

Escolha \(\omega_c\) e \(PM\), e calcule \(\nu_c = \frac{2}{T}\tan(\omega_c T/2)\).
Leia \(\angle G_w(j\nu_c)\) e calcule \(\text{boost} = PM - 90^\circ - \angle G_w(j\nu_c)\), que deve estar em \((0^\circ, 90^\circ)\).
Calcule \(k = \tan\!\left(\dfrac{\text{boost}}{2} + 45^\circ\right)\) e posicione \(\omega_z = \nu_c/k\) e \(\omega_p = k\,\nu_c\).
Ajuste \(K_c\) para que \(|C_w G_w|_{j\nu_c} = 1\).
Obtenha \(C(z) = C_w\big|_{w = \frac{2}{T}\frac{z-1}{z+1}}\).

Função `tipo2_w.m`

function [Cz, info] = tipo2_w(Gz, wc, PM)
T  = Gz.Ts;
Gw = d2c(Gz, 'tustin');                 % plano w
nu = 2/T*tan(wc*T/2);                   % frequência pré-distorcida
[mag, ph] = bode(Gw, nu);               % = G(e^{j wc T})
ph = mod(ph + 180, 360) - 180;
boost = PM - 90 - ph;
assert(boost > 0 && boost < 90, 'Boost fora de (0,90): rever wc/PM')
k  = tand(boost/2 + 45);
wz = nu/k;   wp = nu*k;
s  = tf('s');
Cn = (1 + s/wz)/(s*(1 + s/wp));
Kc = 1/(mag*abs(freqresp(Cn, nu)));
Cz = c2d(Kc*Cn, T, 'tustin');
info = struct('nu',nu,'fase',ph,'boost',boost,'k',k, ...
              'wz',wz,'wp',wp,'Kc',Kc);
end

\(C(z)\) resulta na forma \(K\,\dfrac{(z+1)(z - z_c)}{(z-1)(z-p_c)}\).

Caso 5 — Tipo 2 na planta 1

[C1t2, i1] = tipo2_w(G1z, 20, 60);      % wc = 20 rad/s, PM = 60°
[S, R] = tfdata(C1t2, 'v');

\(\nu_c\)	\(\angle G\)	boost	\(k\)	\(\omega_z\)	\(\omega_p\)	\(K_c\)
20,07	\(-69{,}3^\circ\)	\(39{,}3^\circ\)	2,11	9,52	42,3	10,62

\[ C(z) = \frac{0{,}2042\,z^2 + 0{,}01855\,z - 0{,}1857}{z^2 - 1{,}6507\,z + 0{,}6507} \]

Resultado: \(PM = 60^\circ\), \(M_s = 1{,}42\), sobressinal de 7,7 % e \(t_s \approx 0{,}19\) s.

Caso 6 — Tipo 2 na planta 2

[C2t2, i2] = tipo2_w(G2z, 8, 60);       % wc = 8 rad/s, PM = 60°

\(\nu_c\)	\(\angle G\)	boost	\(k\)	\(\omega_z\)	\(\omega_p\)	\(K_c\)
8,00	\(-68{,}1^\circ\)	\(38{,}1^\circ\)	2,05	3,90	16,4	3,42

Resultado: \(PM = 60^\circ\) ✔, mas \(M_s = 2{,}29\) ✘ (margem de módulo \(0{,}44 < 0{,}5\)), com \(t_s \approx 1{,}9\) s.

Aviso

A margem de fase não basta. Com o cruzamento próximo da ressonância (\(\omega_n = 10\)), a margem de ganho encolhe e \(|S|\) tem um pico alto. Há duas saídas: reduzir \(\omega_c\), o que torna a malha mais lenta, ou acrescentar um segundo zero, passando para o Tipo 3 / PID. Esse é o dilema clássico do buck em modo tensão.

Robustez e comparação

Função `analise.m`

function r = analise(Gz, S, R, nome)
T  = Gz.Ts;
C  = tf(S, R, T);   L = C*Gz;
Sy = feedback(1, L);                          % sensibilidade
F1 = feedback(L, 1);                          % 1 GDL
F2 = minreal(tf(polyval(S,1), R, T)*feedback(Gz, C));  % 2 GDL
[~, Pm, ~, wc] = margin(L);
Ms = norm(Sy, inf);
i1 = stepinfo(F1);   i2 = stepinfo(F2);
r  = table(string(nome), Pm, wc, Ms, i1.Overshoot, i1.SettlingTime, ...
           i2.Overshoot, i2.SettlingTime, 'VariableNames', ...
           {'caso','PM','wc','Ms','OS1','ts1','OS2','ts2'});
figure('Name', nome)
subplot(1,2,1); step(F1, F2); legend('1 GDL','2 GDL'); grid on
subplot(1,2,2); bodemag(Sy); yline(20*log10(2),'--'); grid on
end

Critério de robustez: \(M_s = \|S\|_\infty \le 2\), ou seja, margem de módulo \(\ge 0{,}5\).

Quadro comparativo

Caso	Planta	Método	\(PM\)	\(M_s\)	SS 1 GDL	\(t_s\) 1 GDL	SS 2 GDL	\(t_s\) 2 GDL
PI	1	Diofantina	\(62^\circ\)	1,23	12,8 %	0,16 s	1,5 %	0,12 s
Tipo 2	1	Plano \(w\)	\(60^\circ\)	1,42	7,7 %	0,19 s	—	—
PI	2	LGR	\(92^\circ\)	1,20	≈ 0 %	2,4 s	—	—
Tipo 2	2	Plano \(w\)	\(60^\circ\)	2,29	3,4 %	1,9 s	—	—
PID	2	Diofantina	\(54^\circ\)	1,38	15,6 %	0,37 s	3,9 %	0,46 s

\(t_s\) pelo critério de 2 %. Pequenas diferenças em relação ao stepinfo são esperadas.

Referência: `pidtune` no modelo discreto

C0pi  = pid(1, 1, 'Ts', T, 'IFormula', 'BackwardEuler');
C0pid = pid(1, 1, 1, 0.01, T, 'IFormula','BackwardEuler', ...
                              'DFormula','BackwardEuler');
Cpi  = pidtune(G2z, C0pi);          % PI de referência
Cpid = pidtune(G2z, C0pid, 27);     % PIDF com wc alvo = 27 rad/s
[S, R] = tfdata(tf(Cpid), 'v');
analise(G2z, S, R, 'pidtune PIDF');

O pidtune serve como linha de base: se o projeto manual não supera essa referência em desempenho ou robustez, algo está errado.

Resumo: qual método para qual controlador

Controlador	Método natural	Por quê
PI (planta de 1ª ordem)	Diofantina	solução fechada e exata
PI (planta de ordem superior)	LGR em \(z\) / plano \(w\)	subdeterminado para alocação
PI Tipo 2	Fator \(k\) no plano \(w\)	é um projeto de \(PM\) em \(\omega_c\)
PID	Diofantina (+ 2 GDL)	4 incógnitas para 4 polos (2ª ordem)

Sempre: verifique \(M_s\), porque a alocação de polos e a \(PM\) sozinhas não garantem robustez. Use a estrutura 2 GDL para separar o seguimento de referência da rejeição de perturbação. Respeite a faixa de validade do ARX (banda do PRBS).

Exercícios propostos

Refaça o Caso 3 com \(z_a = 0\) (deadbeat auxiliar). O que acontece com \(p\), \(M_s\) e \(|u|\)?
Projete o Tipo 2 da planta 2 com \(\omega_c = 5\) rad/s. Compare \(M_s\) e \(t_s\) com o Caso 6.
Acrescente um atraso computacional (\(z^{-1}\)) às duas plantas e repita todos os casos. Qual método “enxerga” o atraso sem modificação?
Identifique \(G_2(z)\) com arx a partir de uma PRBS com ruído (SNR = 20 dB) e avalie a degradação de \(M_s\) do PID projetado.

Projeto direto de controladores digitais

Roteiro

Motivação

Emulação × projeto direto

Do ARX ao \(G(z)\)

Plantas e especificações

Plantas de teste (\(T = 10\) ms)

Da especificação ao plano \(z\)

Graus de liberdade: quem resolve o quê

Método 1 — Alocação polinomial

A equação diofantina

Função rst_poly.m

Caso 1 — PI na planta 1 (solução analítica)

Caso 2 — PID na planta 1: o que a álgebra diz

Caso 3 — PID na planta 2

Conversão para ganhos paralelos

Zeros de \(S\) e a estrutura 2 GDL

Implementação no DSP (forma RST)

Método 2 — Lugar das raízes em \(z\)

Caso 4 — PI na planta 2

Leitura do resultado

Método 3 — Resposta em frequência no plano \(w\)

A transformação \(w\)

PI Tipo 2 e o fator \(k\) (Venable)

Função tipo2_w.m

Caso 5 — Tipo 2 na planta 1

Caso 6 — Tipo 2 na planta 2

Robustez e comparação

Função analise.m

Quadro comparativo

Referência: pidtune no modelo discreto

Resumo: qual método para qual controlador

Exercícios propostos

Função `rst_poly.m`

Função `tipo2_w.m`

Função `analise.m`

Referência: `pidtune` no modelo discreto