% Universidade Federal do Rio Grande do Sul
% Escola de Engenharia
% Departamento de Engenharia Elétrica
% ENG04006 - Sistemas e Sinais
%
% Aluno: Alisson Claudino de Jesus (00246796)
% Aluno: Bernardo Brandão Pandolfo (00230340)
% % Turma: C
%
% Laboratório 5 - Questão 1
%
% Considere um sistema representado pela seguinte equação diferencial:
% y'''+41*y''+360*y'+900*y=600*x'+1200*x
% onde y(0)=2, y'(0)=1 e y''(0)=-0.05. Baseado nesse sistema, pede-se:
%
% (a) Calcular a solução dessa equação diferencial, ou seja, calcular a
% resposta natural, a resposta forçada a um degrau unitário e a resposta
% completa, como visto na área 1 da disciplina.
% (b) Determinar a transformada de Laplace do sistema (realizar os cálculos
% à mão ou usando a função 'laplace').
% (c) Através da transformada de Laplace inversa (função 'residue' para
% frações parciais), obter as respostas natural, forçada e completa do
% sistema. Simular e comparar com as respostas obtidas no item (a).
% (d) Usando a transformada de Laplace, simular as resposta natural,
% forçada e completa usando os comandos 'step', 'impulse' ou 'lsim'.
% Comparar com as respostas obtidas anteriormente.
% (e) Implementar um diagrama de blocos que representa a equação
% diferencial dada no simulink (usando blocos somadores, integradores e
% ganhos) e simule as respostas natural, forçada e completa. Obtenha as
% mesmas respostas usando blocos de funções de transferência.
%
% Conclusão: analisando os gráficos correspondentes a cada etapa deste
% trabalho, temos resultados consistentes com aquilo que se é esperado.
%%

clear all;                                                                 %limpa todas as variáveis
close all;                                                                 %fecha todas as janelas

t=0:0.01:10;                                                               %cria o vetor tempo

%(a):
yn1=(0.118.*exp(-30.*t)-13.956.*exp(-6.*t)+15.838.*exp(-5.*t)).*(t>=0);    %resposta natural calculada
yf1=(0.933.*exp(-30.*t)-16.667.*exp(-6.*t)+14.4.*exp(-5.*t)+4/3).*(t>=0);  %resposta forçada calculada
yc1=yn1+yf1;                                                               %resposta completa calculada

figure;                                                                    %abre a janela gráfica 1
subplot(5,1,1);                                                            %formata a janela gráfica com 3 linhas e 1 colunas
plot(t,yn1,'r','linewidth',2);                                             %plota yn em função de t
hold on;                                                                   %mantém o gráfico anterior afim de plotar os gráficos seguintes
plot(t,yf1,'g','linewidth',2);                                             %plota yf em função de t
plot(t,yc1,'b','linewidth',2);                                             %plota yc em função de t
hold off;
axis([0 2 -0.2 3.5]);                                                      %determina os limites do gráfico
ylabel('y(t)');                                                            %nomeia o eixo y
title('Solução calculada');                                                %nomeia a janela gráfica
legend('Resposta natural','Resposta forçada','Resposta completa');         %nomeia os gráficos
grid on;

%(b):
% Y=((2*s^2+83*s+760.95)+(600*s+1200)*X)/(s^3+41*s^2+360*s+900)
% cálculo feito à mão

%(c):
% Resposta natural (x(t)=0, cond. iniciais não-nulas):
% Y=(2*s^2+83*s+760.95)/(s^3+41*s^2+360*s+900)
yn_lap_num=[2 83 760.95];                                                  %coeficientes do numerador de Yn(s)
yn_lap_den=[1 41 360 900];                                                 %coeficientes do denominador de Yn(s)
[an,pn]=residue(yn_lap_num,yn_lap_den);                                    %obtém os coeficientes 'an' e os polos 'pn' das frações parciais

% Resposta forçada (entrada x(t) não-nula, x(0)=0 e cond. iniciais nulas):
% Y=(600*s+1200)/s*(s^3+41*s^2+360*s+900)
yf_lap_num=[600 1200];                                                     %coeficientes do numerador de Yf(s)
yf_lap_den=[1 41 360 900 0];                                               %coeficientes do denominador de Yf(s)
[af,pf]=residue(yf_lap_num,yf_lap_den);                                    %obtém os coeficientes 'af' e os polos 'pf' das frações parciais

% Utilizando os coeficientes encontrados, montamos as respostas natural,
% forçada e completa:
yn2=(an(1).*exp(pn(1).*t)+an(2).*exp(pn(2).*t)+an(3).*exp(pn(3).*t)).*(t>=0);
yf2=(af(1).*exp(pf(1).*t)+af(2).*exp(pf(2).*t)+af(3).*exp(pf(3).*t)+af(4)).*(t>=0);
yc2=yn2+yf2;

subplot(5,1,2);                                                            %formata a janela gráfica com 3 linhas e 1 colunas
plot(t,yn2,'r','linewidth',2);                                             %plota yn em função de t
hold on;                                                                   %mantém o gráfico anterior afim de plotar os gráficos seguintes
plot(t,yf2,'g','linewidth',2);                                             %plota yf em função de t
plot(t,yc2,'b','linewidth',2);                                             %plota yc em função de t
hold off;
axis([0 2 -0.2 3.5]);                                                      %determina os limites do gráfico
ylabel('y(t)');                                                            %nomeia o eixo y
title('Laplace inversa');                                                  %nomeia a janela gráfica
legend('Resposta natural','Resposta forçada','Resposta completa');         %nomeia os gráficos
grid on;

%(d):
% Resposta ao impulso:
num_imp=[2 83 760.95];                                                     %coeficientes do numerador
den_imp=[1 41 360 900];                                                    %coeficientes do denominador
tf1=tf(num_imp,den_imp);                                                   %cria uma função de transferência utilizando os coeficientes listados
yn3=impulse(tf1,t);                                                        %resposta ao impulso da função de transferência criada

% Resposta ao degrau:
num_deg=[600 1200];                                                        %coeficientes do numerador
den_deg=[1 41 360 900];                                                    %coeficientes do denomindador
tf2=tf(num_deg,den_deg);                                                   %cria uma função de transferência utilizando os coeficientes listados
yf3=step(tf2,t);                                                           %resposta ao degrau da função de transferência criada

yc3=yn3+yf3;                                                               %resposta completa utilizando-se os comandos step e impulse

subplot(5,1,3);                                                            %formata a janela gráfica com 3 linhas e 1 colunas
plot(t,yn3,'r','linewidth',2);                                             %plota yn em função de t
hold on;                                                                   %mantém o gráfico anterior afim de plotar os gráficos seguintes
plot(t,yf3,'g','linewidth',2);                                             %plota yf em função de t
plot(t,yc3,'b','linewidth',2);                                             %plota yc em função de t
hold off;
axis([0 2 -0.2 3.5]);                                                      %determina os limites do gráfico
ylabel('y(t)');                                                            %nomeia o eixo y
title('Comandos step e impulse');                                          %nomeia a janela gráfica
legend('Resposta natural','Resposta forçada','Resposta completa');         %nomeia os gráficos
grid on;

%(e):
% Diagramas de blocos:
sim('simulink_q1');

subplot(5,1,4);
plot(tout,yns(:,2),'r','linewidth',2);                                     %tout é o vetor tempo do arquivo simulink
hold on;
plot(tout,yfs(:,2),'g','linewidth',2);
plot(tout,ycs(:,2),'b','linewidth',2);
hold off;
axis([0 2 -0.2 3.5]);                                                      %determina os limites do gráfico
ylabel('y(t)');                                                            %nomeia o eixo y
title('Diagrama de blocos');                                               %nomeia a janela gráfica
legend('Resposta natural','Resposta forçada','Resposta completa');         %nomeia os gráficos
grid on;

% Blocos de função de transferência:
subplot(5,1,5);
plot(tout,yns2(:,2),'r','linewidth',2);
hold on;
plot(tout,yfs2(:,2),'g','linewidth',2);
plot(tout,ycs2(:,2),'b','linewidth',2);
hold off;
axis([0 2 -0.2 3.5]);                                                      %determina os limites do gráfico
ylabel('y(t)');                                                            %nomeia o eixo y
title('Blocos de função de transferência');                                 %nomeia a janela gráfica
legend('Resposta natural','Resposta forçada','Resposta completa');         %nomeia os gráficos
grid on;