%Exercício Proposto 1 - Modelo de planta com incertezas em dois parâmetros.
%Obtenção da função peso para incerteza multiplicativa
%Inicialização
%------------------------------------------------------------------%
clear all
close all
%------------------------------------------------------------------%
%Variáveis do sistema
w=logspace(-2,2); %vetor com a faixa de frequências
am=2; %valor médio de a
bm=4; %valor médio de b
ra=1; % a = am + ra*DELTAa
rb=2;
passo=0.25; % magnitude do "passo" para a variação de DELTA
DELTAa=[-1:passo:1];
DELTAb=[-1:passo:1];
N=size(DELTAa);
N=N(1,2);
M=size(DELTAb);
M=M(1,2);
Gnominal=tf([1],[1 am bm]); %Função de transferência nominal
%------------------------------------------------------------------%
%LOOP responsável pela geração da família de plantas considerada
%O número de curvas geradas será M*N
for i = 1:N
for j = 1:M
DENreal = [1 (am + ra*DELTAa(i)) (bm + rb*DELTAb(j))];
Preal=tf([1],DENreal);
Perro = (Preal-Gnominal)/Gnominal;
[magnitude, fase] = bode(Perro,w);
mg(1:50)=magnitude(1,1,(1:50));
figure(1)
loglog(w,mg);
hold on
end
end
%------------------------------------------------------------------%
%Parte responsável por plotar a curva de W(s)
mf=ginput(5); %Pega 5 pontos no grafico
magg=vpck(mf(:,2),mf(:,1)); %Organiza os pontos para permitir o uso da função "fitmag"
figure(2)
Wim=fitmag(magg); %Escolha da ordem da funçao (CURVE FITTING)
[A,B,C,D]=unpck(Wim) %Converte para espaço de estado

[NUMw,DENw]=ss2tf(A,B,C,D); %Transformação Espaço de estados para FT
W=tf(NUMw,DENw) %FT da função peso (Wim) desejada
save Wsalvo W Gnominal w
%%
%------------------------------------------------------------------%
%Resp em frequência da função peso
load Wsalvo
[magnitudew, fasew] = bode(W,w);
mgw(1:50)=magnitudew(1,1,(1:50));
figure(1)
loglog(w,mgw,'r');

%------------------------------------------------------------------%
%Resp em frequência do sistema com função peso definida
Pmax = Gnominal*(1 + W);
Pmin = Gnominal*(1 - W);
Pint = Gnominal*(1-W*0.6);

[magnitudewmax, fasewmax] = bode(Pmax,w);
mgwmax(1:50)=magnitudewmax(1,1,(1:50));
figure(3)
loglog(w,mgwmax,'r');
hold on

[magnitudewmin, fasewmin] = bode(Pmin,w);
mgwmin(1:50)=magnitudewmin(1,1,(1:50));
figure(3)
loglog(w,mgwmin,'b');
hold on


[magnitudewmed, fasewmed] = bode(Gnominal,w);
mgwmed(1:50)=magnitudewmed(1,1,(1:50));
figure(3)
loglog(w,mgwmed,'g');
hold on

[magnitudeW, faseW] = bode(Pint,w);
mgW(1:50)=magnitudeW(1,1,(1:50));
figure(3)
loglog(w,mgW,'c');
legend('Parcela Máx', 'Parcela Mín', 'Nominal', 'Intermediária');
%%