%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Projeto  de  Auxílio  à  Pesquisa  Fo-
% mentado   Pela   Fundação   de   Amparo
% à  Pesquisa  do  Estado  de  São  Paulo
% (FAPESP), processo 2019/14317-2

% Desenvolvimento de um software para estimação de
% parâmetros de máquinas elétricas pelos ensaios de
% rejeição de carga.

%Bolsista: Pedro Frizzo de Moura
%Orientador: Elmer Pablo Tito Cari


%Última alteração no software: 26/02/2020
%Criador: Pedro Frizzo de Moura
%Colaborador: Bruno Fioravante


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Algoritmo para encontrar parâmetros de eixo direto para o teste de rejeição de carga

close all
clear all
clc

%% ------------------------------------------------------------------------------- %
% Máquina do artigo  (BORTONI; JARDINI, 2002) 
% 9375 kVA / 13.8 kV / 60 Hz 
Uo=0.972; % Valor de pico da tensão de armadura antes do curto
Xd = 1.1375;      % Xd
Xid = 0.1876;     % X'd
Xiid = 0.1074;    % X''d
Tid = 4.9653;     % T'd
Tiid = 0.0222;    % T''d
Xq = 1.055;
Xiiq=0.1492;
Tiiq=0.074;
beta = 0;
f = 60;         % Hz
t=0:0.0001:40;
w=2*pi*f;
E = 1; %pu
iq0 = (Uo * sin(beta))/Xq;
id0 = (E - (Uo * cos(beta))) / Xd;
id0 = 0.321;
for i=1:length(t)
     Ua(i) = Uo * sin(w*t(i) - beta)  - iq0 *[Xq - (Xq - Xiiq)*exp(-t(i)/Tiiq)]*cos(w*t(i)) - id0*[Xd - (Xd - Xid)*exp(-t(i)/Tid)-(Xid - Xiid)*exp(-t(i)/Tiid)]*sin(w*t(i));
end

 %Função para encontra máximos e mínimos
[picos,posicoes] = findpeaks(Ua,t,'MinPeakDistance',0.003); 
[vales,locais]=findlows(Ua,t);

figure(1)
plot(t,Ua)
hold on
plot(posicoes,picos,'o')
hold on
plot(locais,vales,'*')
legend('Ua(t)','máximos','mínimos'); xlabel('tempo (s)'); ylabel('tensão (pu)');
title('Ua(t) e seus máximos e mínimos');

%% Obtendo a média das envoltórias
[x,y]=media_env(picos,vales,posicoes,locais);

figure(2)
plot(x,y,'.')
ylabel('tensão(pu)'); xlabel('tempo(s)'); title('Média das envoltórias --> udec(t)')
var = [(y(1) - y(30)) / y(30)] * 20;
[valor_regime,instante,int]=regime(y,x,var)
Vregime=valor_regime;

Xd_found=(Uo-Vregime)/id0;

% Retirando o valor de regime permanente da média da envoltória
figure(3)
v2 = y-Vregime;
plot(x,v2)
ylabel('tensão(pu)'); xlabel('tempo(s)'); title('udec(t)-Uregime')

% Retirando o regime do vetor
v3=v2(1:int); % Sinal com os regimes transitorio e subtransitorio
x3=x(1:int);    % Vetor com tempos de 0 até "instante que atinge o regime"

figure(4)
plot(x3,v3)  
ylabel('tensão(V)'); xlabel('tempo(s)'); title('Apenas as partes transitória e subtransitória do sinal "udec(t) - Uregime"')

v3_log=log(v3);

figure(5)
plot(x3,v3_log,'.')
ylabel('log[udec(t)-Uregime]'); xlabel('tempo(s)'); title('"udec(t) - Uregime" em escala logarítmica')
hold on
%% 
% PARA ACHAR O PONTO ONDE TERMINA O REGIME SUBTRANSIENTE NA ESCALA LOG

x4=x3(length(x3)/10:end);
v4=v3_log(length(v3)/10:end);
%%
% Encontrando coeficientes da reta mais aproximada
coef=polyfit(x4,v4,1);
a1=coef(2);
b1=coef(1);

% Gerando a reta
Reta_trans = a1 + b1*x3;
% Plotando a reta em cima dos pontos de log
plot(x3,Reta_trans)

Xid_found= Xd - (exp(a1)/id0);
Tid_found= -1/b1;

%% Retirando também o regime transitório de "udec(t)-Vregime" que vai de t=0 até o tempo que entra em regime
% Primeiro, gera-se os valores de regime transitório
valor_trans= id0*(Xd_found - Xid_found) * exp(-x3/Tid_found);

% Retirando o transitório, ponto a ponto, de "um(t) - Vregime"
v5=v3-valor_trans;

figure(7)
plot(x3,v5, 'o');
% plot(x3,v3)  
% hold on
% plot(x3,v5,'*')
% legend('um(t)-Vregime','um(t)-Vregime-transitório')
% Observa-se que o regime subtransitorio corresponde aos 4 primeiros dados
% do vetor, futuramente o reconhecimento dos regimes será feito de forma
% automática através do refinamento de código.

v5_log=log(v5);

x6=x3(1:4);
v6=v5_log(1:4);
figure(8)
plot(x3, v5_log, 'o');
% Encontrando coeficientes da reta mais aproximada
coef2=polyfit(x6,v6,1);
a2=coef2(2);
b2=coef2(1);

% Gerando a reta
Reta_subtrans = a2 + b2*x3;

figure(9)
plot(x6,v6,'o');
hold on;
plot(Reta_subtrans,'--');

Xiid_found= Xid_found - (exp(a2)/id0);
Tiid_found= -1/b2;


