matlab_EV充电功率计算 电动汽车充电功率

% Newton Raphson LoadFlow % Adilson Carlos Batista % adilsoncbatista@gmail.com % https://www.facebook.com/dilsinho.ee clear all;close all;clc tic% contando o tempo. te = @() datestr(now);% hora e data da simula玢o. disp(te()) ;% informando hora e data. format short c = 0 ;% contador para o n鷐ero de itera玢o. c1 = 1 ;% contador para auxiliar a encontrar barra em paralelo %% Par鈓etros da Linha de Transmiss鉶 %ev-plug in电动汽车接入潮流计算 t=18; n=1500; nev=zeros(1,24); for i=1:24 nev(i)=n*(normpdf(i,17.6,3.4)+normpdf(i+24,17.6,3.4)); end nEV=zeros(24,10000); for i=1:24 nEV(i,:)=poissrnd(nev(i),[1 10000])*3.6/1000; end nu=zeros(1,24); no=zeros(1,24); for i=1:24 uev(i)=mean(nEV(i,:)); oev(i)=var(nEV(i,:)).^0.5; evsk(i)=skewness(nEV(i,:)); evku(i)=kurtosis(nEV(i,:)); wev1(i)=evsk(i)*(1/2)+(evku(i)-(3/4)*(evsk(i))^2)^0.5; wev2(i)=evsk(i)*(1/2)-(evku(i)-(3/4)*(evsk(i))^2)^0.5; wev3=0; end Pev=zeros(24,3); for i=1:24 Pev(i,1)=uev(i)+oev(i)*wev1(i); Pev(i,2)=uev(i)+oev(i)*wev2(i); Pev(i,3)=uev(i)+oev(i)*wev3; end wwe(1)=(1/1)*(wev1(t)*(wev1(t)-wev2(t)))^(-1); wwe(2)=-(1/1)*(wev2(t)*(wev1(t)-wev2(t)))^(-1); wwe(3)=1/1-wwe(1)-wwe(2); V=zeros(3,33); for i=1:3 % load('dados_da_linha_ex1_ASEP1.mat','LineData') % LineData = [de para R X Bsh tap] LineData = [ 1 2 0.012 0.033 0 1 2 3 0.078 0.211 0 1 3 4 0.048 0.130 0 1 4 5 0.05 0.136 0 1 5 6 0.055 0.149 0 1 6 7 0.064 0.0173 0 1 7 8 0.12 0.324 0 1 8 9 0.116 0.314 0 1 9 10 0.124 0.336 0 1 10 11 0.105 0.283 0 1 11 12 0.102 0.276 0 1 12 13 0.101 0.273 0 1 13 14 0.071 0.193 0 1 14 15 0.078 0.210 0 1 15 16 0.098 0.266 0 1 16 17 0.091 0.245 0 1 17 18 0.083 0.225 0 1 2 19 0.127 0.343 0 1 19 20 0.198 0.535 0 1 20 21 0.12 0.323 0 1 21 22 0.106 0.288 0 1 3 23 0.138 0.374 0 1 23 24 0.171 0.462 0 1 24 25 0.118 0.319 0 1 6 26 0.027 0.072 0 1 26 27 0.037 0.101 0 1 27 28 0.014 0.038 0 1 28 29 0.093 0.251 0 1 29 30 0.067 0.181 0 1 30 31 0.062 0.169 0 1 31 32 0.041 0.112 0 1 32 33 0.071 0.192 0 1]; [nLine,nCol] = size(LineData) ;% N鷐ero de linhas do sistema. %% Dados das Barras do Sistema El閠rico de Pot阯cia Sb = 100 ; % Pot阯cia base para todo sistema (omitindo MVA para todas pot阯cias). Vb = 230 ; % Tens鉶 base para todo sistema (omitindo KV para todas tens鮡s). Zb = Vb^2 / Sb ; % Imped鈔cia base para todo sistema (Ohms) % bd = [barra tipo Pgi Qgi Pli Qli Qgmin Qgmax Pgmin Pgmax Vesp Vmin Vmax Angulo] ; BusData = [ 1 1 0 0 0 0 0 0 0 0 1.05 0.9 1.05 0; 2 3 0 0 0.266 0.178 -10 10 10 50 1.000 .8 1.05 0 ; 3 3 0 0 0.200 0.133 -40 50 30 250 1.000 .8 1.05 0 ; 4 3 0 0 0.244 0.178 -10 10 10 50 1.000 .8 1.05 0 ; 5 3 0 0 0.222 0.155 0 0 0 0 1.000 .8 1.06 0 ; 6 3 0 0 0.266 0.200 0 0 0 0 1.000 .8 1.05 0 ; 7 3 0 0 0.155 0.111 -6 24 0 0 1.000 .8 1.05 0 ; 8 3 0 0 0.333 0.222 0 0 0 0 1.000 .8 1.05 0 ; 9 3 0 0 0.289 0.178 -6 24 0 0 1.000 .8 1.05 0 ; 10 3 0 0 0.222 0.155 0 0 0 0 1.000 .8 1.05 0 ; 11 3 0 0 0.222 0.155 0 0 0 0 1.000 .8 1.05 0 ; 12 3 0 0 0.211 0.155 0 0 0 0 1.000 .8 1.06 0 ; 13 3 0 0 0.444 0.333 0 0 0 0 1.000 .8 1.06 0 ; 14 3 0 0 0.400 0.333 0 0 0 0 1.000 .8 1.06 0 ; 15 3 0 0 0.289 0.222 0 0 0 0 1.000 .8 1.06 0 ; 16 3 0 0 0.333 0.266 0 0 0 0 1.000 .8 1.05 0 ; 17 3 0 0 0.333+Pev(t,i) 0.222 0 0 0 0 1.000 .8 1.06 0 ; 18 3 0 0 0.222 0.178 0 0 0 0 1.000 .8 1.06 0 ; 19 3 0 0 0.289 0.178 0 0 0 0 1.000 .8 1.06 0 ; 20 3 0 0 0.333 0.266 0 0 0 0 1.000 .8 1.06 0 ; 21 3 0 0 0.333 0.222 0 0 0 0 1.000 .8 1.05 0 ; 22 3 0 0 0.444 0.133 0 0 0 0 1.000 .8 1.06 0 ; 23 3 0 0 0.444 0.222 0 0 0 0 1.000 .8 1.06 0 ; 24 3 0 0 0.333 0.266 0 0 0 0 1.000 .8 1.06 0 ; 25 3 0 0 0.266 0.222 0 0 0 0 1.000 .8 1.06 0 ; 26 3 0 0 0.178 0.133 0 0 0 0 1.000 .8 1.05 0 ; 27 3 0 0 0.222 0.111 0 0 0 0 1.000 .8 1.06 0 ; 28 3 0 0 0.222 0.155 0 0 0 0 1.000 .8 1.06 0 ; 29 3 0 0 0.266 0.155 0 0 0 0 1.000 .8 1.06 0 ; 30 3 0 0 0.289 0.155 0 0 0 0 1.000 .8 1.06 0 ; 31 3 0 0 0.155 0.089 0 0 0 0 1.0 0.8 1.05 0; 32 3 0 0 0.155 0.089 -10 10 10 50 1.000 .8 1.05 0 ; 33 3 0 0 0.133 0.089 -40 50 30 250 1.000 .8 1.05 0 ]; bus = BusData(:,1); % Numeros que Identificam as Barras nbus = max(bus) ; % N鷐ero de baras do sistema. type = BusData(:,2) ; % Tipo para as barras 1-Vteta, 2-PV, 3-PQ Pg = BusData(:,3)/Sb ; % PGi Qg = BusData(:,4)/Sb ; % QGi Pl = BusData(:,5)/Sb ; % PLi Ql = BusData(:,6)/Sb ; % QLi Qgmin = BusData(:,7)/Sb ; % Limite mininimo de pot阯cia reativa Qgmax = BusData(:,8)/Sb ; % Limite maximo de pot阯cia reativa Pgmin = BusData(:,9)/Sb ; % Limite mininimo de pot阯cia ativa Pgmax = BusData(:,10)/Sb ; % Limite maximo de pot阯cia ativa Vesp = BusData(:,11) ; % Tens鉶 inicial especificada. Vmin = BusData(:,12) ; % Tens鉶 Minima Permitida na Barra. Vmax = BusData(:,13) ; % Tens鉶 Maxima Permitida na Barra. Ang = BusData(:,14) ; % Angulo inicial especificado f1 = find(type == 1) ; % Identificando a barra de refer阯cia f2 = find(type == 2) ; % Identificando as barras do tipo PV. f3 = find(type == 3) ; % Identificando as barras do tipo PQ. busSW = f1' ; % Barra de refer阯cia. busPV = f2' ; % Vetor das Barras PV . busPQ = f3' ; % Vetor das Barras PQ . NPV = length(busPV) ; % Numero de Barras PV. NPQ = length(busPQ) ; % Numero de Barras PQ. busPVPQ = sort([busPQ busPV]) ; % Indices dos angulos theta a serem encontrados. %% Prealoca玢o de Mem髍ia Ybus = zeros(nbus,nbus) ;% Matriz das Admit鈔cia nas Barras. H = zeros(nbus,nbus) ;% Acoplamento da pot阯cia ativa com o 鈔gulo theta. M = zeros(nbus,nbus) ;% Acoplamento da pot阯cia reativa com o 鈔gulo theta. N = zeros(nbus,nbus) ;% Acoplamento da pot阯cia ativa com a tens鉶 V. L = zeros(nbus,nbus) ;% Acoplamento da pot阯cia reativa com a tens鉶 V. Pcal = zeros(nbus,nbus) ;% Pot阯cia ativa a ser calculada.(forma Matricial) Qcal = zeros(nbus,nbus) ;% Pot阯cia ativa a ser calculada.(forma Matricial) Pf = zeros(2*nLine,1) ;% Fluxo de pot阯cia ativa. Qf = zeros(2*nLine,1) ;% Fluxo de pot阯cia reativa. Ppl = zeros(nLin