%% Projeto de um Controlador de Avanço-Atraso de Fase em MATLAB
% Limpa o workspace e a janela de comando
clear;
clc;
% Define 's' como uma variável de função de transferência
s = tf('s');
%% 1. Definição do Processo P(s)
% P(s) = 192.1 / (s * (s+6) * (s+10))
P_num = 192.1;
P_den = [1 16 60 0]; % Coeficientes de s^3 + 16s^2 + 60s
P = tf(P_num, P_den);
disp('----------------------------------------------');
disp('1. Processo P(s):');
disp(P);
%% 2. Especificações de Projeto
zeta_d = 0.456; % Fator de amortecimento desejado
wn_d = 7.869; % Frequência natural não amortecida desejada
ess_ramp_d = 0.0313; % Erro em regime permanente para rampa unitária desejado
disp('----------------------------------------------');
disp('2. Especificações de Projeto:');
fprintf(' Fator de Amortecimento (zeta) = %.3f\n', zeta_d);
fprintf(' Frequência Natural (wn) = %.3f rad/s\n', wn_d);
fprintf(' Erro em Regime Permanente (ess para rampa) = %.4f\n', ess_ramp_d);
%% 3. Projeto do Compensador de Avanço (Lead Compensator)
% 3.1. Cálculo dos Polos Dominantes Desejados (sd)
% sd = -zeta*wn +/- j*wn*sqrt(1-zeta^2)
sd_real = -zeta_d * wn_d;
sd_imag = wn_d * sqrt(1 - zeta_d^2);
sd = sd_real + 1i * sd_imag; % Polo superior no plano complexo
disp('----------------------------------------------');
disp('3. Projeto do Compensador de Avanço:');
fprintf(' Polos Desejados (sd) = %.4f +/- j%.4f\n', sd_real, sd_imag);
% 3.2. Parâmetros do Compensador de Avanço (calculados previamente)
Kc = 10.295;
zc = 6; % Zero do compensador de avanço (para cancelar polo em s=-6)
pc = 29.1555; % Polo do compensador de avanço
C_lead = Kc * (s + zc) / (s + pc);
disp(' Compensador de Avanço C_lead(s):');
disp(C_lead);
%% 4. Projeto do Compensador de Atraso (Lag Compensator)
% 4.1. Cálculo do Kv Requerido
Kv_required = 1 / ess_ramp_d;
disp('----------------------------------------------');
disp('4. Projeto do Compensador de Atraso:');
fprintf(' Kv Requerido = %.4f\n', Kv_required);
% 4.2. Cálculo do Kv do Sistema com Compensador de Avanço (G_OL_lead)
% G_OL_lead = C_lead * P;
% Para calcular Kv, precisamos da função de malha aberta C_lead * P
% Observação: O polo em s=-6 de P(s) é cancelado pelo zero de C_lead(s)
G_OL_lead_simplified_num = Kc * P_num; % 10.295 * 192.1 = 1977.8
G_OL_lead_simplified_den = [1 10] .* [1 29.1555] .* [1 0]; % s * (s+10) *
(s+29.1555)
G_OL_lead_simplified_den = conv([1 10], [1 29.1555]); % (s+10)(s+29.1555)
G_OL_lead_simplified_den = conv(G_OL_lead_simplified_den, [1 0]); % s * (s+10)
(s+29.1555)
% G_OL_lead_simplified = tf(G_OL_lead_simplified_num, G_OL_lead_simplified_den);
% disp(G_OL_lead_simplified); % Mostra a função de transferência simplificada
% A maneira robusta de calcular Kv é usando o limite de s*G_OL(s) quando s->0.
% Ou de forma mais prática, obter o ganho DC de s*G_OL_lead, após eliminar polos em
s=0.
% Para um sistema Tipo 1, Kv = lim_{s->0} s*G_OL(s).
% Então, Kv = (Kc * P_num * z_l) / (p_l * 10 * 29.1555) se o polo em 0 não fosse
cancelado.
% A forma exata é: Kv = dcgain(s * C_lead * P) se s*C_lead*P não tiver polo em 0.
% Usando a expressão analítica do Kv para o sistema com C_lead (após cancelamento
de s+6):
% Kv_lead = (Kc * P_num) / (10 * 29.1555);
% Kv_lead = (10.295 * 192.1) / (10 * 29.1555);
% Kv_lead = 1977.8 / 291.555;
% Kv_lead = 6.7844; % Usando o valor exato da etapa de projeto
% Em MATLAB, para um sistema Tipo 1, K_v = lim_{s->0} s * G(s)
% Se G(s) = N(s) / (s * D(s)), então s * G(s) = N(s) / D(s).
% Kv = dcgain(N(s)/D(s))
G_OL_lead_unsimplified = C_lead * P;
[num_ol_lead, den_ol_lead] = tfdata(G_OL_lead_unsimplified, 'v');
% Encontra o polo em s=0 no denominador e o remove simbolicamente
% Exemplo: se den = [1 10 0], remover o 0 para obter [1 10]
Kv_lead_calc = num_ol_lead(end) / (den_ol_lead(end-1) * den_ol_lead(end));
% Esta é uma forma simplificada, mas a mais robusta é usar dcgain(s*G_OL)
Kv_lead_matlab = dcgain(s * G_OL_lead_unsimplified);
fprintf(' Kv do sistema com apenas Compensador de Avanço = %.4f\n',
Kv_lead_matlab);
% Nota: Pequenas diferenças de arredondamento podem ocorrer entre o cálculo manual
e MATLAB.
% 4.3. Parâmetros do Compensador de Atraso (calculados previamente)
zl = 0.1; % Zero do compensador de atraso
pl = 0.0212; % Polo do compensador de atraso
C_lag = (s + zl) / (s + pl);
disp(' Compensador de Atraso C_lag(s):');
disp(C_lag);
%% 5. Controlador Final C(s)
C = C_lead * C_lag;
disp('----------------------------------------------');
disp('5. Controlador Final C(s) = C_lead(s) * C_lag(s):');
disp(C);
%% 6. Verificação do Projeto
% 6.1. Função de Malha Aberta Total
G_OL_total = C * P;
% 6.2. Função de Malha Fechada (assumindo realimentação unitária)
G_CL_total = feedback(G_OL_total, 1);
disp('----------------------------------------------');
disp('6. Verificação do Projeto:');
% 6.3. Verificação dos Polos de Malha Fechada e Lugar das Raízes
figure;
rlocus(G_OL_total);
hold on;
plot(real(sd), imag(sd), 'rx', 'MarkerSize', 10, 'LineWidth', 2, 'DisplayName',
'Polos Desejados');
plot(real(sd), -imag(sd), 'rx', 'MarkerSize', 10, 'LineWidth', 2);
sgrid(zeta_d, wn_d); % Desenha linhas de amortecimento e frequência natural
title('Lugar das Raízes do Sistema Compensado');
legend('show', 'Location', 'best');
grid on;
hold off;
cl_poles = pole(G_CL_total);
disp(' Polos de Malha Fechada do Sistema Compensado:');
disp(cl_poles);
% Espera-se que haja polos próximos a 'sd' e polos adicionais introduzidos pelo
compensador.
% Os polos próximos à origem (do lag compensator) e o polo de alta frequência (do
lead compensator)
% são tipicamente não-dominantes.
% 6.4. Verificação da Resposta Transitória (Resposta ao Degrau)
figure;
step(G_CL_total);
title('Resposta ao Degrau do Sistema de Malha Fechada Compensado');
grid on;
% Obtém informações detalhadas da resposta ao degrau
info = stepinfo(G_CL_total);
disp(' Informações da Resposta ao Degrau (para verificação de Zeta e Wn):');
fprintf(' Overshoot: %.2f%%\n', [Link]);
fprintf(' Tempo de Assentamento (2%%): %.2f s\n', [Link]);
fprintf(' Pico: %.4f\n', [Link]);
fprintf(' Tempo de Pico: %.2f s\n', [Link]);
% Note: stepinfo calcula essas métricas com base na resposta global.
% A validação de zeta e wn é primariamente pelo posicionamento dos polos.
% 6.5. Verificação do Erro em Regime Permanente para Rampa
% Kv = lim_{s->0} s * G_OL(s)
% Calculando o ganho DC de s * G_OL_total(s)
Kv_final = dcgain(s * G_OL_total);
ess_ramp_final = 1 / Kv_final;
disp(' Verificação do Erro em Regime Permanente para Rampa:');
fprintf(' Constante de Velocidade (Kv) Final = %.4f\n', Kv_final);
fprintf(' Erro em Regime Permanente (ess) Final para Rampa = %.4f\n',
ess_ramp_final);
fprintf(' Erro Requerido (ess para rampa) = %.4f\n', ess_ramp_d);
disp('----------------------------------------------');
disp('Projeto do Controlador Concluído!');
disp('Verifique os gráficos e os valores para confirmar que as especificações foram
atendidas.');
disp('Os polos de malha fechada devem estar próximos aos polos desejados.');
disp('O Kv final deve ser próximo ao Kv requerido, resultando no ess desejado.');