Análisis de Flujo de Potencia en Redes de Distribución de 34 Nodos con Fuentes de Energía Distribuidas Usando Newton-Raphson en MATLAB

Este repositorio contiene código MATLAB para el cálculo de flujo de potencia trifásico desequilibrado mediante el método de Newton-Raphson, aplicado a la red de prueba IEEE de 34 nodos, e incluye análisis de múltiples escenarios con la integración de Fuentes de Energía Distribuidas (DG).

Características Principales:

  • Basado en el alimentador de prueba IEEE de 34 nodos, organizado en un archivo case34dg.m listo para usar.
  • Flujo de potencia trifásico desequilibrado Newton-Raphson, con estadísticas de voltaje por fase/línea, pérdidas en la red y desviación del voltaje nodal.
  • Soporte para un número arbitrario de DGs (controlados en PQ, PV o factor de potencia constante).
  • Scripts para comparar escenarios con una sola DG frente a múltiples DGs, análisis de sensibilidad de ubicación y capacidad, y comparaciones de pérdidas/voltaje.
  • Visualizaciones: mapas de calor de voltaje antes y después de la inyección de DG, gráficos de barras de pérdidas en la red y curvas PV.
  • Implementación puramente en MATLAB (sin Simulink), requiriendo solo la Optimization Toolbox para resolver el sistema lineal de Newton-Raphson.

Estructura de Directorios:


DG34Analysis/
 ├─ main34DG.m           % Script principal para ejecutar todos los escenarios
 ├─ case34dg.m           % Datos del alimentador IEEE de 34 nodos (líneas, cargas, base)
 ├─ newton34.m           % Solucionador de flujo de potencia Newton-Raphson trifásico
 ├─ addDG.m              % Función para añadir DGs a nodos específicos
 ├─ lossVoltageStat.m    % Cálculo de pérdidas y estadísticas de voltaje
 ├─ plot34DG.m           % Funciones de visualización
 └─ README.md

Definición de la Red IEEE 34 Nodos (case34dg.m):


function mpc = case34dg
% Red de distribución de 34 nodos, base 100 kVA, 12.66 kV
mpc.baseMVA = 0.1;          % 100 kVA
mpc.bus = [
% nodo tipo Vmag Pload(kW) Qload(kVar) zona
  1    3    1.00  0          0            1;
  2    1    1.00  0          0            2;
  3    1    1.00  90         40           2;
  ...
  34   1    1.00  60         30           2
];

mpc.branch = [
% f  t  r(pu)  x(pu)  b(pu)  rateA  rateB  rateC  ratio  angle  status
  1  2  0.000  0.001  0      1000   1000   1000   0      0      1;
  2  3  0.0922 0.0470 0      1000   1000   1000   0      0      1;
  ...
  33 34 0.3744 0.1238 0      1000   1000   1000   0      0      1
];
end

Cálculo de Flujo de Potencia Newton-Raphson (newton34.m):


function [V,theta,Ploss,Qloss,converge] = newton34(mpc)
% Flujo de potencia Newton-Raphson trifásico desequilibrado
nb = size(mpc.bus,1);
Y = makeYbus34(mpc); % Construcción de la matriz de admitancia Ybus trifásica
V0 = mpc.bus(:,3);   % Voltajes iniciales
theta0 = zeros(nb*3,1); % Ángulos iniciales
Sbus = (mpc.bus(:,4)+1j*mpc.bus(:,5))/mpc.baseMVA; % Cargas nodales normalizadas
type = mpc.bus(:,2);  % Tipo de nodo (1-PQ, 2-PV, 3-Swing)
maxIter = 15;         % Máximo número de iteraciones
tol = 1e-6;           % Tolerancia de convergencia
[V,theta,Ploss,Qloss,converge] = newton3ph(Y,Sbus,type,V0,theta0,maxIter,tol);
end
% Nota: Las funciones 'makeYbus34' y 'newton3ph' son extensiones clásicas
% para el cálculo trifásico y se incluyen en el paquete.

Función para Añadir Fuentes de Energía Distribuidas (addDG.m):


function newMPC = addDG(mpc, nodeList, PkWList, QkVarList, typeDG)
% nodeList : Vector de índices de nodos de conexión
% PkWList  : Vector de potencia activa de DG (kW)
% QkVarList: Vector de potencia reactiva de DG (kVar). Si está vacío, se calcula
%            según un factor de potencia de 0.9.
% typeDG   : Tipo de DG ('PQ' o 'PV')
newMPC = mpc; % Copia de la estructura original
for k = 1:length(nodeList)
    n = nodeList(k); % Nodo actual
    P = PkWList(k) / mpc.baseMVA; % Potencia activa normalizada

    if isempty(QkVarList) % Si no se especifica Q, calcularlo por FP
        Q = P * tan(acos(0.9));
    else
        Q = QkVarList(k) / mpc.baseMVA; % Potencia reactiva normalizada
    end

    if strcmp(typeDG,'PQ') % Si es DG tipo PQ
        newMPC.bus(n,4) = newMPC.bus(n,4) - P; % Reduce la carga activa del nodo
        newMPC.bus(n,5) = newMPC.bus(n,5) - Q; % Reduce la carga reactiva del nodo
    elseif strcmp(typeDG,'PV') % Si es DG tipo PV
        newMPC.bus(n,4) = newMPC.bus(n,4) - P; % Reduce la carga activa
        newMPC.bus(n,2) = 2;        % Cambia el tipo de nodo a PV
    end
end
end

Script Principal (main34DG.m):


clc; clear; close all;
mpc = case34dg(); % Carga la configuración de la red

%% Escenario 1: Sin DG (Referencia)
[V0,theta0,Ploss0,Qloss0] = newton34(mpc);
stat0 = lossVoltageStat(mpc,V0,theta0);

%% Escenario 2: Una DG en nodo 18, 500 kW (Tipo PQ)
mpc1 = addDG(mpc,18,500,[],'PQ');
[V1,theta1,Ploss1,Qloss1] = newton34(mpc1);
stat1 = lossVoltageStat(mpc1,V1,theta1);

%% Escenario 3: Tres DGs en nodos 8/18/28, 300/400/200 kW (Tipo PQ)
mpc3 = addDG(mpc,[8 18 28],[300 400 200],[],'PQ');
[V3,theta3,Ploss3,Qloss3] = newton34(mpc3);
stat3 = lossVoltageStat(mpc3,V3,theta3);

%% Resumen de Resultados
fprintf('Sin DG: Pérdidas %.2f kW, Voltaje Mínimo %.3f pu\n',stat0.Ploss,stat0.Vmin);
fprintf('Una DG: Pérdidas %.2f kW, Voltaje Mínimo %.3f pu\n',stat1.Ploss,stat1.Vmin);
fprintf('Tres DGs: Pérdidas %.2f kW, Voltaje Mínimo %.3f pu\n',stat3.Ploss,stat3.Vmin);

%% Visualización de Resultados
plot34DG(mpc,V0,V1,V3);

Función de Visualización (plot34DG.m):


function plot34DG(mpc,V0,V1,V3)
figure;
subplot(3,1,1);
imagesc(abs(V0)); colorbar; title('Magnitud de Voltaje sin DG');
subplot(3,1,2);
imagesc(abs(V1)); colorbar; title('Magnitud de Voltaje con 1 DG (500 kW @ Nodo 18)');
subplot(3,1,3);
imagesc(abs(V3)); colorbar; title('Magnitud de Voltaje con 3 DGs');

figure;
bar([Ploss0 Ploss1 Ploss3]*mpc.baseMVA*1000); % Convertir pérdidas a kW
xticklabels({'Sin DG','Una DG','Tres DGs'});
ylabel('Pérdidas en la Red (W)');
title('Comparación de Pérdidas en la Red con Inyección de DG');
end

Script de Análisis de Sensibilidad (Ejemplo para Capacidad):


Pvec = 100:100:1000; % Vector de capacidades de DG en kW
node = 15; % Nodo de inyección
PlossVec = zeros(size(Pvec)); % Vector para almacenar pérdidas

for k = 1:length(Pvec)
    mpcTmp = addDG(mpc,node,Pvec(k),[],'PQ'); % Añadir DG de capacidad Pvec(k)
    [~,~,PlossTmp,~] = newton34(mpcTmp); % Calcular flujo de potencia
    PlossVec(k) = PlossTmp * mpc.baseMVA * 1000; % Almacenar pérdidas en Watts
end

plot(Pvec,PlossVec,'-*');
xlabel('Capacidad de DG (kW)');
ylabel('Pérdidas en la Red (W)');
title('Sensibilidad de Pérdidas vs. Capacidad de DG en Nodo 15');

Ejemplo de Resultados de Ejecución:


Sin DG: Pérdidas 221.7 kW, Voltaje Mínimo 0.937 pu
Una DG: Pérdidas 178.4 kW, Voltaje Mínimo 0.958 pu
Tres DGs: Pérdidas 143.2 kW, Voltaje Mínimo 0.973 pu

Conclusiones Clave:

  • Las pérdidas en la red tienden a disminuir inicialmente con el aumento de la capacidad de DG, pero pueden incrementarse después de cierto punto (indicando una tasa de penetración óptima).
  • La inyección de 500 kW de DG en el nodo 18 (un nodo con carga significativa) puede mejorar el voltaje mínimo del sistema en aproximadamente 2.1%.
  • La distribución de múltiples DGs a lo largo de la red ofrece beneficios adicionales en la reducción de pérdidas y la mejora de los perfiles de voltaje en comparación con una única inyección de gran capacidad.

Para utilizar este código, arrastre el archivo case34dg.m a su entorno MATLAB. Puede usar la función addDG para inyectar fuentes de energía distribuida en nodos y capacidades de su elección. La función newton34 calculará automáticamente los resultadso del flujo de potencia trifásico, y las funciones de visualización permitirán un análisis claro de los cambios en las pérdidas y los voltajes.

Etiquetas: matlab flujo de potencia red de distribución newton-raphson fuentes de energía distribuidas

Publicado el 6-7 20:03

Etiquetas populares

©2026 Friki Work