% ATV03 - Fluxo de Potencia
% Aluno: Lucas Monte Bravo
% Programa: GNU Octave / MATLAB compativel
%
% IMPORTANTE NO OCTAVE ONLINE:
% 1) Cole este codigo inteiro em um arquivo chamado, por exemplo:
%       Fluxo_Potencia_Lucas_Monte_Bravo.m
% 2) Clique em Run.
% 3) Nao execute por partes no console.

clear;
clc;
format long;

fprintf('ATV03 - Fluxo de Potencia\n');
fprintf('Aluno: Lucas Monte Bravo\n');
fprintf('Programa: GNU Octave / MATLAB compativel\n\n');

% ============================================================
% DADOS DO ALUNO
% ============================================================
P2 = 0.04;          % Potencia ativa gerada na barra 2, PV
P3c = 0.10;         % Carga ativa na barra 3
Q3c = 0.02;         % Carga reativa na barra 3
L = 0.03;

% Convencao usada no fluxo de potencia:
% injecao positiva; carga negativa.
P3 = -P3c;
Q3 = -Q3c;

V1 = 1.00;
theta1 = 0.00;

V2esp = 0.98;       % Barra 2 PV

% ============================================================
% MATRIZ YBUS
% ============================================================
Y = zeros(3,3);

% Dados das linhas: [de para r x bsh]
% bsh ja e metade da susceptancia da linha, conforme enunciado.
linhas = [
    1 2 2*L 20*L 0.2*L;
    1 3 4*L 40*L 0.4*L;
    2 3 2*L 20*L 0.2*L
];

for k = 1:size(linhas,1)
    i = linhas(k,1);
    j = linhas(k,2);
    r = linhas(k,3);
    x = linhas(k,4);
    bsh = linhas(k,5);

    y = 1/(r + 1i*x);

    Y(i,i) = Y(i,i) + y + 1i*bsh;
    Y(j,j) = Y(j,j) + y + 1i*bsh;

    Y(i,j) = Y(i,j) - y;
    Y(j,i) = Y(j,i) - y;
end

fprintf('Matriz Ybus:\n');
disp(Y);

% ============================================================
% 1) GAUSS-SEIDEL - 3 ITERACOES
% ============================================================
fprintf('\n============================================================\n');
fprintf('1) GAUSS-SEIDEL - 3 ITERACOES\n');
fprintf('============================================================\n');

% Flat start:
% Barra 1: slack, 1 ang 0
% Barra 2: PV, magnitude especificada 0.98 ang 0
% Barra 3: PQ, 1 ang 0
V = [1 + 0i; V2esp + 0i; 1 + 0i];

fprintf('Iteracao        theta2(graus)        V3(pu)        theta3(graus)        |dV2|        |dV3|\n');

maiorErro3 = 0;

for it = 1:3
    Vant = V;

    % --------------------------
    % Barra 2 - PV
    % --------------------------
    Icalc = Y*V;
    Scalc = V .* conj(Icalc);
    Q2calc = imag(Scalc(2));

    soma = Y(2,1)*V(1) + Y(2,3)*V(3);

    V2novo = ((P2 - 1i*Q2calc)/conj(V(2)) - soma)/Y(2,2);

    % Mantem |V2| = 0.98 pu
    V(2) = V2esp * V2novo/abs(V2novo);

    % --------------------------
    % Barra 3 - PQ
    % --------------------------
    soma = Y(3,1)*V(1) + Y(3,2)*V(2);

    V(3) = ((P3 - 1i*Q3)/conj(V(3)) - soma)/Y(3,3);

    dV2 = abs(V(2) - Vant(2));
    dV3 = abs(V(3) - Vant(3));

    theta2GS = angle(V(2))*180/pi;
    V3GS = abs(V(3));
    theta3GS = angle(V(3))*180/pi;

    fprintf('%4d          %14.6f      %10.6f       %14.6f      %8.6f    %8.6f\n', ...
        it, theta2GS, V3GS, theta3GS, dV2, dV3);

    if it == 3
        maiorErro3 = max(dV2,dV3);
    end
end

fprintf('\nMaior erro da terceira iteracao = %.6f pu\n', maiorErro3);

% ============================================================
% 2) NEWTON-RAPHSON
% ============================================================
fprintf('\n============================================================\n');
fprintf('2) NEWTON-RAPHSON\n');
fprintf('============================================================\n');

% Variaveis de estado:
% x = [theta2; theta3; V3]
theta2 = 0;
theta3 = 0;
V3 = 1;

tol = 0.003;
maxIter = 20;

% minIter = 2 deixa a tabela com duas linhas.
% Se quiser parar exatamente na primeira vez que epsilon < 0.003,
% altere para minIter = 1.
minIter = 2;

fprintf('Iteracao        theta2(graus)        V3(pu)        theta3(graus)        P1(pu)        Q1(pu)        Q2(pu)        epsilon\n');

for it = 1:maxIter

    % --------------------------------------------------------
    % Calcula P e Q no estado atual
    % --------------------------------------------------------
    Vmag = [V1; V2esp; V3];
    theta = [theta1; theta2; theta3];

    Pcalc = zeros(3,1);
    Qcalc = zeros(3,1);

    for i = 1:3
        for k = 1:3
            Gik = real(Y(i,k));
            Bik = imag(Y(i,k));
            ang = theta(i) - theta(k);

            Pcalc(i) = Pcalc(i) + Vmag(i)*Vmag(k)*(Gik*cos(ang) + Bik*sin(ang));
            Qcalc(i) = Qcalc(i) + Vmag(i)*Vmag(k)*(Gik*sin(ang) - Bik*cos(ang));
        end
    end

    % Mismatch:
    % Barra 2 PV: DeltaP2
    % Barra 3 PQ: DeltaP3 e DeltaQ3
    mismatch = [
        P2 - Pcalc(2);
        P3 - Pcalc(3);
        Q3 - Qcalc(3)
    ];

    % --------------------------------------------------------
    % Jacobiano
    % Linhas: DeltaP2, DeltaP3, DeltaQ3
    % Colunas: theta2, theta3, V3
    % --------------------------------------------------------
    J = zeros(3,3);

    % Linha 1: P2, barra i = 2
    i = 2;

    % dP2/dtheta2
    J(1,1) = -Qcalc(i) - imag(Y(i,i))*Vmag(i)^2;

    % dP2/dtheta3
    k = 3;
    Gik = real(Y(i,k));
    Bik = imag(Y(i,k));
    ang = theta(i) - theta(k);
    J(1,2) = Vmag(i)*Vmag(k)*(Gik*sin(ang) - Bik*cos(ang));

    % dP2/dV3
    J(1,3) = Vmag(i)*(Gik*cos(ang) + Bik*sin(ang));

    % Linha 2: P3, barra i = 3
    i = 3;

    % dP3/dtheta2
    k = 2;
    Gik = real(Y(i,k));
    Bik = imag(Y(i,k));
    ang = theta(i) - theta(k);
    J(2,1) = Vmag(i)*Vmag(k)*(Gik*sin(ang) - Bik*cos(ang));

    % dP3/dtheta3
    J(2,2) = -Qcalc(i) - imag(Y(i,i))*Vmag(i)^2;

    % dP3/dV3
    J(2,3) = Pcalc(i)/Vmag(i) + real(Y(i,i))*Vmag(i);

    % Linha 3: Q3, barra i = 3
    i = 3;

    % dQ3/dtheta2
    k = 2;
    Gik = real(Y(i,k));
    Bik = imag(Y(i,k));
    ang = theta(i) - theta(k);
    J(3,1) = -Vmag(i)*Vmag(k)*(Gik*cos(ang) + Bik*sin(ang));

    % dQ3/dtheta3
    J(3,2) = Pcalc(i) - real(Y(i,i))*Vmag(i)^2;

    % dQ3/dV3
    J(3,3) = Qcalc(i)/Vmag(i) - imag(Y(i,i))*Vmag(i);

    % --------------------------------------------------------
    % Resolve J*dx = mismatch
    % --------------------------------------------------------
    dx = J\mismatch;

    theta2 = theta2 + dx(1);
    theta3 = theta3 + dx(2);
    V3 = V3 + dx(3);

    % --------------------------------------------------------
    % Recalcula P e Q apos a atualizacao
    % --------------------------------------------------------
    Vmag = [V1; V2esp; V3];
    theta = [theta1; theta2; theta3];

    Pcalc = zeros(3,1);
    Qcalc = zeros(3,1);

    for i = 1:3
        for k = 1:3
            Gik = real(Y(i,k));
            Bik = imag(Y(i,k));
            ang = theta(i) - theta(k);

            Pcalc(i) = Pcalc(i) + Vmag(i)*Vmag(k)*(Gik*cos(ang) + Bik*sin(ang));
            Qcalc(i) = Qcalc(i) + Vmag(i)*Vmag(k)*(Gik*sin(ang) - Bik*cos(ang));
        end
    end

    epsilon = max(abs([
        P2 - Pcalc(2);
        P3 - Pcalc(3);
        Q3 - Qcalc(3)
    ]));

    fprintf('%4d          %14.6f      %10.6f       %14.6f      %10.6f    %10.6f    %10.6f    %10.8f\n', ...
        it, theta2*180/pi, V3, theta3*180/pi, Pcalc(1), Qcalc(1), Qcalc(2), epsilon);

    if epsilon < tol && it >= minIter
        break;
    end
end

fprintf('\nResultado final usado:\n');
fprintf('theta2 = %.6f graus\n', theta2*180/pi);
fprintf('V3     = %.6f pu\n', V3);
fprintf('theta3 = %.6f graus\n', theta3*180/pi);
fprintf('P1     = %.6f pu\n', Pcalc(1));
fprintf('Q1     = %.6f pu\n', Qcalc(1));
fprintf('Q2     = %.6f pu\n', Qcalc(2));
fprintf('epsilon final = %.8f\n', epsilon);

fprintf('\nTabela final das barras:\n');
fprintf('Barra    Tipo       P(pu)        Q(pu)        V(pu)        theta(graus)\n');
fprintf('  1      Slack   %10.6f  %10.6f  %10.6f      %10.6f\n', Pcalc(1), Qcalc(1), V1, 0);
fprintf('  2      PV      %10.6f  %10.6f  %10.6f      %10.6f\n', P2, Qcalc(2), V2esp, theta2*180/pi);
fprintf('  3      PQ      %10.6f  %10.6f  %10.6f      %10.6f\n', P3, Q3, V3, theta3*180/pi);