%% Resolução do Problema: Projeto de Controlador PID em MATLAB
% Limpa o workspace e o command window
clear;
clc;
%% 1. Definição dos Polos Dominantes Desejados
% Especificações de desempenho
Tp = 1; % Tempo de pico (segundos)
zeta = 0.7; % Relação de amortecimento
% Cálculo dos parâmetros do polo dominante
wd = pi / Tp; % Frequência natural amortecida
wn = wd / sqrt(1 - zeta^2); % Frequência natural não amortecida
sigma = zeta * wn; % Componente real do polo
% Polos dominantes desejados
s_d = -sigma + 1i*wd; % Um dos polos (o outro é o conjugado)
fprintf('1. Polos Dominantes Desejados:\n');
fprintf(' sigma = %.4f\n', sigma);
fprintf(' wd = %.4f rad/s\n', wd);
fprintf(' wn = %.4f rad/s\n', wn);
fprintf(' s_d = %.4f + j%.4f\n\n', real(s_d), imag(s_d));
%% 2. Análise da Planta Original
% Função de transferência da planta G(s) = K / ((s+1.2)(s+3.7))
% Vamos assumir K=1 para o projeto do controlador, e o ganho K_D do PID
% irá absorver o ganho total necessário para atender o critério de magnitude.
K_plant = 1;
num_G = K_plant;
den_G = conv([1, 1.2], [1, 3.7]); % Multiplica (s+1.2) por (s+3.7)
G_s = tf(num_G, den_G);
fprintf('2. Planta Original G(s):\n');
% CORREÇÃO APLICADA AQUI (se já não estiver): Usando disp()
disp(G_s); % Use disp() para exibir a função de transferência
fprintf('\n');
%% 3. Projeto do Controlador PID usando o Critério do Lugar das Raízes
% Polos do sistema em malha aberta (incluindo o polo do integrador do PID)
p_plant = pole(G_s); % Obtém os polos da planta diretamente
p_integrator = 0; % Polo do integrador do PID
all_poles = [p_plant; p_integrator];
% Calcular os ângulos dos polos em relação ao polo desejado s_d
angles_poles = angle(s_d - all_poles) * (180/pi); % Convertendo para graus
fprintf('3. Projeto do Controlador PID:\n');
fprintf(' a) Critério de Ângulo para os Zeros do PID:\n');
fprintf(' Ângulos dos polos em relação a s_d (em graus):\n');
fprintf(' p1 (-1.2): %.2f\n', angles_poles(1));
fprintf(' p2 (-3.7): %.2f\n', angles_poles(2));
fprintf(' p3 (0): %.2f\n', angles_poles(3));
sum_angles_poles = sum(angles_poles);
fprintf(' Soma dos ângulos dos polos: %.2f graus\n', sum_angles_poles);
% O ângulo que os zeros do PID devem fornecer
% Sum(angles_zeros) - Sum(angles_poles) = (2k+1)*180
% Para k=0, Sum(angles_zeros) = 180 + Sum(angles_poles)
required_angle_zeros = 180 + sum_angles_poles;
fprintf(' Ângulo total requerido dos zeros do PID (base): %.2f degrees\n',
required_angle_zeros);
% Normalizando o ângulo para o intervalo principal (-180, 180]
% Ajuste para garantir que o ângulo esteja na faixa esperada para tan().
while required_angle_zeros > 180
required_angle_zeros = required_angle_zeros - 360;
end
while required_angle_zeros <= -180
required_angle_zeros = required_angle_zeros + 360;
end
fprintf(' Ângulo efetivo dos zeros (para cálculo): %.2f degrees\n',
required_angle_zeros);
% Se os dois zeros do PID são iguais (-z0_positive), então cada zero contribui com
metade do ângulo
angle_per_zero = required_angle_zeros / 2;
fprintf(' Ângulo por zero do PID: %.2f degrees\n', angle_per_zero);
% Calcular a posição do zero (-z0_positive)
% z0_positive = sigma + wd / tan(angle_per_zero_rad)
% Atenção: o cálculo original manual era z0 = 3,079 + 0,7126 = 3,7916
% O zero está em -z0_positive. Então, a distância real de s_d até o zero é
(z0_positive - sigma)
% Ou seja, (real(s_d) - (-z0_positive)) = (z0_positive - sigma)
% tan(angle_per_zero_rad) = imag(s_d) / (z0_positive - sigma)
% z0_positive - sigma = imag(s_d) / tan(angle_per_zero_rad)
% z0_positive = sigma + imag(s_d) / tan(angle_per_zero_rad)
z0_positive_calculated = -real(s_d) + imag(s_d) / tand(angle_per_zero); % A forma
mais robusta
fprintf(' Posição calculada para o zero duplo (-z0_positive): %.4f\n', -
z0_positive_calculated);
% b) Critério de Magnitude para o Ganho K_D
% |C(s_d) * G(s_d)| = 1
% C(s) = K_D * (s + z0_positive_calculated)^2 / s
% G(s) = K_plant / ((s + 1.2)(s + 3.7))
% Módulos dos termos no ponto s_d
abs_sd = abs(s_d);
abs_sd_plus_1_2 = abs(s_d + 1.2);
abs_sd_plus_3_7 = abs(s_d + 3.7);
abs_sd_plus_z0_sq = abs(s_d + z0_positive_calculated)^2;
K_D = (abs_sd * abs_sd_plus_1_2 * abs_sd_plus_3_7) / (abs_sd_plus_z0_sq * K_plant);
fprintf('\n b) Critério de Magnitude para o Ganho K_D:\n');
fprintf(' Módulos no ponto s_d:\n');
fprintf(' |s_d| = %.4f\n', abs_sd);
fprintf(' |s_d+1.2| = %.4f\n', abs_sd_plus_1_2);
fprintf(' |s_d+3,7| = %.4f\n', abs_sd_plus_3_7);
fprintf(' |s_d + %.4f|^2 = %.4f\n', z0_positive_calculated,
abs_sd_plus_z0_sq);
fprintf(' Ganho K_D do PID = %.4f\n', K_D);
% Determinar Kp, Ki, Kd a partir de K_D e z0_positive_calculated
% C(s) = K_D * (s^2 + 2*z0_positive_calculated*s + z0_positive_calculated^2) / s
% C(s) = K_D*s + (2*z0_positive_calculated*K_D) + (z0_positive_calculated^2*K_D)/s
Kd_final = K_D;
Kp_final = 2 * z0_positive_calculated * K_D;
Ki_final = z0_positive_calculated^2 * K_D;
fprintf('\n Ganhos do Controlador PID:\n');
fprintf(' Kp = %.4f\n', Kp_final);
fprintf(' Ki = %.4f\n', Ki_final);
fprintf(' Kd = %.4f\n', Kd_final);
%% 4. Controlador PID Proposto
% Criar o objeto PID no MATLAB
% Num: [Kd Kp Ki] para Kds^2 + Kps + Ki
% Den: [1 0] para s
C_s_pid_coeffs = tf([Kd_final, Kp_final, Ki_final], [1, 0]);
% Ou na forma de zero-polo-ganho (zpk)
zeros_C = [-z0_positive_calculated, -z0_positive_calculated]; % Zeros são -
z0_positive_calculated
poles_C = 0; % Polo em s=0
C_s_zpk = zpk(zeros_C, poles_C, Kd_final); % Ganho K_D
(mesmo que Kd_final)
fprintf('\n4. Controlador PID Proposto:\n');
% CORREÇÃO APLICADA AQUI: Usando disp() para exibir os objetos tf e zpk
fprintf(' C(s) (coeficientes):\n');
disp(C_s_pid_coeffs);
fprintf(' C(s) (na forma ZPK):\n');
disp(C_s_zpk);
fprintf('\n');
%% 5. Função de Transferência de Malha Fechada
% L(s) = C(s) * G(s)
L_s = C_s_pid_coeffs * G_s; % Usamos a forma de coeficientes para a multiplicação
% T(s) = L(s) / (1 + L(s)) para realimentação unitária negativa
T_s = feedback(L_s, 1);
fprintf('\n5. Função de Transferência de Malha Aberta L(s):\n');
disp(L_s); % Usando disp()
fprintf('\n5. Função de Transferência de Malha Fechada T(s):\n');
disp(T_s); % Usando disp()
fprintf('\n');
%% Verificação da Resposta ao Degrau e Lugar das Raízes
fprintf('\nVerificando a Resposta ao Degrau:\n');
figure;
step(T_s);
grid on;
title('Resposta ao Degrau do Sistema em Malha Fechada com Controlador PID');
xlabel('Tempo (segundos)');
ylabel('Amplitude');
% Obter informações da resposta ao degrau para verificar as especificações
step_info = stepinfo(T_s);
fprintf(' Tempo de Pico (Tp): %.4f segundos\n', step_info.PeakTime);
fprintf(' Sobressinal (Overshoot): %.2f%%\n', step_info.Overshoot);
% Calcular zeta a partir do overshoot percentual (se houver overshoot)
if step_info.Overshoot > 0
MO_ratio = step_info.Overshoot / 100;
% Evitar log(0) se MO_ratio for muito pequeno (quase sem overshoot)
if MO_ratio > 1e-6 % Um pequeno limiar para evitar problemas numéricos
ln_MO_ratio = log(MO_ratio);
zeta_from_MO = sqrt( (ln_MO_ratio)^2 / (pi^2 + (ln_MO_ratio)^2) );
fprintf(' Zeta estimado do Overshoot: %.4f\n', zeta_from_MO);
else
fprintf(' Sobressinal muito baixo, zeta próximo de 1 ou maior.\n');
end
else
fprintf(' Não houve sobressinal (Overshoot <= 0).\n');
end
% Plot do Lugar das Raízes com os polos desejados
figure;
rlocus(L_s); % Plota o lugar das raízes da função de malha aberta L(s)
hold on;
plot(real(s_d), imag(s_d), 'rx', 'MarkerSize', 10, 'LineWidth', 2, 'DisplayName',
'Polo Desejado'); % Polo desejado
plot(real(s_d), -imag(s_d), 'rx', 'MarkerSize', 10, 'LineWidth', 2,
'HandleVisibility', 'off'); % Conjugado (não repetir na legenda)
% Adiciona os polos e zeros de malha aberta para melhor visualização
plot(real(all_poles), imag(all_poles), 'bo', 'MarkerSize', 8, 'LineWidth', 1.5,
'DisplayName', 'Polos de Malha Aberta');
plot(real(zeros_C), imag(zeros_C), 'gs', 'MarkerSize', 8, 'LineWidth', 1.5,
'DisplayName', 'Zeros do Controlador');
title('Lugar das Raízes com Controlador PID e Polos Desejados');
xlabel('Eixo Real');
ylabel('Eixo Imaginário');
% Ajusta a legenda para evitar duplicatas e ser mais clara
legend('show', 'Location', 'best');
grid on;
axis equal; % Garante que os eixos tenham a mesma escala para uma visualização
correta do ângulo
sgrid(zeta, wn); % Plota as linhas de damping ratio e frequência natural
hold off;