Диссертация (1144079), страница 13
Текст из файла (страница 13)
5clear; clc;format long e;EPS = [0 0.005 0.01 0.025];for ie = 1:length(EPS)eps = EPS(ie);fname = sprintf('Err_One%d', 1000*eps);load(fname,'-mat');clear fname;Fz = 18; %Set FontSizefigure;ax = gca;grid on;surfl(x1,x2,abs(I1));xlabel('x1 [Ohm]');ylabel('x2 [Ohm]');zlabel('I1 [A]');%ax.View = [65 30];ax.FontSize = Fz;fname = sprintf('Ieps%d', 1000*eps);99%print(fname, '-dpng', '-r300')clear fname;figure;ax = gca;grid on;surfl(x1,x2,100*abs(abs(I1)-abs(Curr1))./abs(Curr1));xlabel('x1 [Ohm]');ylabel('x2 [Ohm]');zlabel('Error of I1 (%)');%ax.View = [65 30];ax.FontSize = Fz;fname = sprintf('ErrI_One%d', 1000*eps);print(fname, '-dpng', '-r300')clear fname ax;figure;ax = gca;grid on;surfl(x1,x2,abs(U1));xlabel('x1 [Ohm]');ylabel('x2 [Ohm]');zlabel('U1 [V]');%ax.View = [65 30];ax.FontSize = Fz;fname = sprintf('Ueps%d', 1000*eps);%print(fname, '-dpng', '-r300')clear fname ax;figure;ax = gca;grid on;surfl(x1,x2,100*abs(abs(U1)-abs(Vol1))./abs(Vol1));xlabel('x1 [Ohm]');ylabel('x2 [Ohm]');zlabel('Error of U1 [%]');%ax.View = [65 30];ax.FontSize = Fz;fname = sprintf('ErrU_One%d', 1000*eps);print(fname, '-dpng', '-r300')clear fname ax;Epsilon_I_everEpsilon_U_everendПрил.
6clear; clc;format long e;EPS = [0 0.005 0.01 0.025];for ie = 1:length(EPS)eps = EPS(ie);fname = sprintf('Aver_Err_Mul%d.mat', 1000*eps);load(fname)100clear fname;Fz = 18; %Set FontSizefigure;ax = gca;grid on;surfl(x1,x2,erI1ever);xlabel('x1 [Ohm]');ylabel('x2 [Ohm]');zlabel('Average error of I1 (%)');%ax.View = [65 30];ax.FontSize = Fz;fname = sprintf('Aver_ErrI_Mul%d', 1000*eps);print(fname, '-dpng', '-r300')clear fname;figure;ax = gca;grid on;surfl(x1,x2,erU1ever);xlabel('x1 [Ohm]');ylabel('x2 [Ohm]');zlabel('Average error of U1 (%)');ax.FontSize = Fz;fname = sprintf('Aver_ErrU_Mul%d', 1000*eps);print(fname, '-dpng', '-r300')clear fname;Epsilon_I_everEpsilon_U_everendПрил.
7clear; clc;format long e;EPS = [0 0.005 0.01 0.025];BUS = [1 2 4 5 6];WeiCoe = [ 1 1 sqrt(2); 1 1 2];XMAX1 = [0.5 1 1 1 1 2];XMAX2 = [0.5 1 1 1 1 2];TIMES = 1;for nBus = 1%1:length(BUS)Bus = BUS(nBus ); % 1 2 4 5 6;for nTimes = 1:TIMESTimes = nTimes;for WeC = 1:2for ie = 1%1:length(EPS)SymSec2i;endendend101endПрил. 8format long e; profile on;tab = fopen('circuit.txt');TAB = fscanf(tab, '%d %d %d %d %d %d %d %d %d %d %d %d', [13 inf]);fclose(tab);A = TAB.';NYG = 1; nBranch = (size(A,2)-NYG)/3; nNode = (size(A,1)-NYG)/3;tab = fopen('local.txt');Lo = fscanf(tab, '%d %d', [2 inf]);fclose(tab);Lo=Lo.';Xmax=[XMAX1(Bus) XMAX2(Bus)];% limit of reactive conductivityXmin=[0 0];eps = EPS(ie);% loadfname = sprintf('RealData2%d.txt', Bus);tab = fopen(fname);clear fname;colon = char(58);[FG,count1] = fscanf(tab,['%d' colon '%d' '%f' '%f' '%f' '%f' '%f' '%f' '%f' '%f''%f' '%f' '%f' '%f' '%f' '%f' '%f' '%f' '%f' '%f'],[20 inf]);fclose(tab);nvalue = 300;for p = 1:720%size(FG,2)ZA(1:nBranch)=[1e-9;0.001063+j*0.005905;0.8*FG(18,p)*(FG(15,p)+j*sin(acos(FG(15,p)))); 1.31429-j*8.89675];ZB(1:nBranch)=[1e-9;0.001063+j*0.005905;0.7*FG(19,p)*(FG(16,p)+j*sin(acos(FG(16,p)))); 1.31429-j*8.89675];ZC(1:nBranch)=[1e-9;0.001063+j*0.005905;FG(20,p)*(FG(17,p)+j*sin(acos(FG(17,p)))); 1.31429-j*8.89675];if NYG == 0ZN = [];elseZN(1:NYG) = 0.0017+j*0.0089;% if the circuit contain Zn!!!endZ=[ZA ZB ZC ZN];Y=1./Z;clear ZA ZB ZC ZN;Y = diag(Y);E = zeros(size(A,2),1);102E(1,1) = 220;%FG(3,p);E(5,1) = 220*exp(-2*pi*i/3);%FG(4,p)*exp(-2*pi*i/3);E(9,1) = 220*exp(2*pi*i/3);%FG(5,p)*exp(2*pi*i/3);ncase = 1;nvar = 2;for ncase = 1:size(Lo,1)X1 = [0 0.1 0.5 1 1.5 0.2 0.4 1.3];X2 = [1.5 1 0.8 0.7 0.4 0.1 0.5 0];for nex = 1:size(X1,2)Z1 = Z;Z1(Lo(ncase,1)) = Z1(Lo(ncase,1)) + j*X1(nex);Z1(Lo(ncase,2)) = Z1(Lo(ncase,2)) + j*X2(nex);Y1 = diag(1./Z1);% Solve circuitU0 = -(A*Y1*A')\(A*Y1*E);Us = A'*U0;U = Us+E;I = Y1*U;Uex(:,nex) = U.*(1 + eps*2*(rand(size(A,2),1) - 0.5));Iex(:,nex) = I.*(1 + eps*2*(rand(size(A,2),1) - 0.5));end% tim he so cua bieu thuc dong dien nhanh chua X1for nf = 1:size(X1,2)Vx(nf,:)=[1X2(nf)-Iex(Lo(ncase,1),nf)*X1(nf)Iex(Lo(ncase,1),nf)*X2(nf) -Iex(Lo(ncase,1),nf)*X1(nf)*X2(nf)];Tx(nf,1) = [Iex(Lo(ncase,1),nf)];end;Ki(:,Lo(ncase,1)) = (Vx.'*Vx)\(Vx.'*Tx);clear Vx Tx nf;% tim he so cua bieu thuc dong dien nhanh chua X2for nf = 1:3Vy(nf,:) = [1 X1(nf)];Ty(nf,1)=[(1+Ki(3,Lo(ncase,1))*X1(nf)+Ki(4,Lo(ncase,1))*X2(nf)+Ki(5,Lo(ncase,1))*X1(nf)*X2(nf))*Iex(Lo(ncase,2),nf)];end;Ki(1:2,Lo(ncase,2)) = (Vy.'*Vy)\(Vy.'*Ty);clear Vy Ty nf;% tim he so cua bieu thuc dong dien cac nhanh con laifor ni = 1:size(A,2)if ni == Lo(ncase,1) || ni == Lo(ncase,2)continueelsefor nf = 1:6103V(nf,:) = [1 X1(nf) X2(nf) X1(nf)*X2(nf)];T(nf,1)=[(1+Ki(3,Lo(ncase,1))*X1(nf)+Ki(4,Lo(ncase,1))*X2(nf)+Ki(5,Lo(ncase,1))*X1(nf)*X2(nf))*Iex(ni,nf)];end;Ki(1:4,ni) = (V.'*V)\(V.'*T);clear nf V T;endend; clear nifor nu = 1:size(A,2)for nf = 1:6V(nf,:) = [1 X1(nf) X2(nf) X1(nf)*X2(nf)];T(nf,1)=[(1+Ki(3,Lo(ncase,1))*X1(nf)+Ki(4,Lo(ncase,1))*X2(nf)+Ki(5,Lo(ncase,1))*X1(nf)*X2(nf))*Uex(nu,nf)];end;Ku(1:4,nu) = (V.'*V)\(V.'*T);clear nf V T;end; clear nu;% Mathematical Optimization% First step, we choose 'n_row' random value of com.dev.
and then wecan get variable power flow, correlatively.n_row=1000;for n = 1:nvarData(:,n)=Xmin(n)+(Xmax(n)-Xmin(n))*rand(n_row,1);endclear n;for n = 1:n_row% define variable power flows% define branch's voltagefor nn = 1:size(A,2)Data(n,nvar+nn)=(Ku(1,nn)+Ku(2,nn)*Data(n,1)+Ku(3,nn)*Data(n,2)+Ku(4,nn)*Data(n,1)*Data(n,2))/(1+Ki(3,Lo(ncase,1))*Data(n,1)+Ki(4,Lo(ncase,1))*Data(n,2)+Ki(5,Lo(ncase,1))*Data(n,1)*Data(n,2));endclear nn;% define branch's currents104for nb=1:size(A,2)if nb == Lo(ncase,1)Data(n,nvar+size(A,2)+nb)=(Ki(1,nb)+Ki(2,nb)*Data(n,2))/(1+Ki(3,Lo(ncase,1))*Data(n,1)+Ki(4,Lo(ncase,1))*Data(n,2)+Ki(5,Lo(ncase,1))*Data(n,1)*Data(n,2));elseif nb == Lo(ncase,2)Data(n,nvar+size(A,2)+nb)=(Ki(1,nb)+Ki(2,nb)*Data(n,1))/(1+Ki(3,Lo(ncase,1))*Data(n,1)+Ki(4,Lo(ncase,1))*Data(n,2)+Ki(5,Lo(ncase,1))*Data(n,1)*Data(n,2));elseData(n,nvar+size(A,2)+nb)=(Ki(1,nb)+Ki(2,nb)*Data(n,1)+Ki(3,nb)*Data(n,2)+Ki(4,nb)*Data(n,1)*Data(n,2))/(1+Ki(3,Lo(ncase,1))*Data(n,1)+Ki(4,Lo(ncase,1))*Data(n,2)+Ki(5,Lo(ncase,1))*Data(n,1)*Data(n,2));endend; clear nbData(n,nvar+2*size(A,2)+1:nvar+3*size(A,2))=Data(n,nvar+1:nvar+size(A,2)).*Data(n,nvar+size(A,2)+1:nvar+2*size(A,2));%PvarData(n,nvar+3*size(A,2)+1)=sum(Data(n,nvar+2*size(A,2)+1:nvar+3*size(A,2)));%sum PvarS = [1 1 1; exp(-2*pi*i/3) exp(2*pi*i/3) 1; exp(2*pi*i/3) exp(2*pi*i/3) 1];Data(n,nvar+3*size(A,2)+2:nvar+3*size(A,2)+4) = (S\[Data(n,16);Data(n,20); Data(n,24)]).';% I secquenceData(n,nvar+3*size(A,2)+5:nvar+3*size(A,2)+7) = (S\[Data(n,17);Data(n,21); Data(n,25)]).';Data(n,nvar+3*size(A,2)+8:nvar+3*size(A,2)+10) = (S\[Data(n,18);Data(n,22); Data(n,26)]).';Data(n,nvar+3*size(A,2)+11:nvar+3*size(A,2)+13) = (S\[Data(n,19);Data(n,23); Data(n,27)]).';endnloop = 8;alpha(1)=1;beta=(10/0.3)^(1/nloop);for k=2:nloopalpha(k)=alpha(k-1)/beta;endgamma = alpha/5;for n = 1:nloop105%po = Data(:,nvar+3*size(A,2)+8:nvar+3*size(A,2)+10);for m = 1:1:n_rowpz1 = 0;for nl = 1:1:n_row% find number of point belong to zone 1if(abs(Data(nl,nvar+3*size(A,2)+10))<abs(Data(m,nvar+3*size(A,2)+10))) & ((1/abs(Data(nl,nvar+3*size(A,2)+8))) <(1/abs(Data(m,nvar+3*size(A,2)+8)))) & (abs(Data(nl,nvar+3*size(A,2)+9)) <abs(Data(m,nvar+3*size(A,2)+9)))pz1 = pz1+1;endend; clear nl;Data(m,nvar+3*size(A,2)+14) = pz1;end; clear m ;% get nvalue valueSpz = 0;INDpz = find(Data(:,nvar+3*size(A,2)+14) <= Spz);while size(INDpz,1) < nvalueSpz = Spz+1;INDpz = find(Data(:,nvar+3*size(A,2)+14) <= Spz);end;if size(INDpz,1) == nvalueData = Data(INDpz,:);elseData = Data(INDpz,:);Data_temp = abs(Data);Data_temp(:,size(Data_temp,2)+1) = sqrt((Data_temp(:,49)).^2+(Data_temp(:,50)).^2+(Data_temp(:,51)).^2);% eliminate redundant valuesif Spz == 0% get 5 rows, which 1/I1 approach minData_temp = sortrows(Data_temp,[-49 51 50]);Data_temp1 = Data_temp(1:5,:);Data_temp(1:5,:) = [];% get 5 rows, which I2 approach minData_temp = sortrows(Data_temp,[50 51 -49]);Data_temp2 = Data_temp(1:5,:);Data_temp(1:5,:) = [];% get 5 rows, which I0 approach minData_temp = sortrows(Data_temp,[51 -49 50]);Data_temp3 = Data_temp(1:5,:);Data_temp(1:5,:) = [];106% get 5 rows which sqrt((1/I1)^2+I2^2+I0^2) approach minData_temp = sortrows(Data_temp,size(Data_temp,2));Data_temp4 = Data_temp(1:5,:);Data_temp(1:5,:) = [];% get ..
more rows in the remaining rows of Data_tempn_part=size(Data_temp,1)/(nvalue-size(Data_temp1,1)size(Data_temp2,1)-size(Data_temp3,1));for n_p = 1:(nvalue-size(Data_temp1,1)-size(Data_temp2,1)size(Data_temp3,1))m_p = fix((n_p-1)*n_part + ceil(n_part*rand));Data_temp5(n_p,:) = Data_temp(m_p,:);end% nvalue rowsData=[Data_temp1(:,1:56);Data_temp2(:,1:56);Data_temp3(:,1:56); Data_temp4(:,1:56); Data_temp5(:,1:56)];elseINDpz_less = find(Data(:,nvar+3*size(A,2)+14) <= (Spz-1));Data_temp1 = Data_temp(INDpz_less,:);Data_temp(INDpz_less,:) = [];Data_temp2 = Data_temp(ceil(size(Data_temp,1)*rand(nvaluesize(INDpz_less,1),1)),:);Data = [Data_temp1(:,1:56); Data_temp2(:,1:56)];endend; clear Data_temp Data_temp1 Data_temp2 Data_temp3Data_temp4 Data_temp5 INDpz_less;Data(:,nvar+3*size(A,2)+15) = 1;%./Data(:,nvar+3*size(A,2)+14);%Data(:,nvar+3*size(A,2)+15)=sqrt(abs(Data(:,nvar+3*size(A,2)+10)).^2+abs(Data(:,nvar+3*size(A,2)+8)).^-2);%domi(:,:,n,p)=Data;divpart=Data(1:nvalue,nvar+3*size(A,2)+15)/sum(Data(1:nvalue,nvar+3*size(A,2)+15));test(1)=0;for k=1:length(divpart)test(k+1)=test(k)+(k);endclear k;h=nvalue;while h<n_rowt=rand;107k=1;while ((t<test(k)) || (t>=test(k+1)))%k=k+1;endt=ceil(nvar*rand);% change x1 x2 randomlytempo = Data(k,:);tempo(t) = tempo(t)*(1+alpha(n)*2*(rand-0.5));if(tempo(t) < Xmin(t)|| tempo(t) > Xmax(t))% limit of x1,x2,...