Help me port this Matlab code......pls....(ported)[Solved]
11-25-2015, 05:16 PM (This post was last modified: 11-29-2015 07:22 AM by toshk.)
Post: #1
 toshk Member Posts: 192 Joined: Feb 2015
Help me port this Matlab code......pls....(ported)[Solved]
Code:
 linedata=[1     2       0.01008,    0.05040,   3.815629,     -19.078144,     10.25,  0.05125;           1     3       0.00744,    0.03720,   5.169561,     -25.847809,     7.75,   0.03875;           2     4       0.00744,    0.03720,   5.169561,     -25.847809,     7.75,   0.03875;           3     4       0.01272,    0.06360,   3.023705,     -15.118528,     12.75,  0.06375];                  busdata=[1   0,     0,  50,     30.99,  1.00,   0   1;          2   0,     0,  170,    105.35, 1.00,   0   2;          3   0,     0,  200,    123.94, 1.00,   0   2;          4   318,   0,   80,    49.58,  1.02,   0   3] % Bus Type: 1.Slack Bus    2.PQ Bus    3.PV Bus       ss=1i*linedata(:,8);      y=linedata(:,5)+1i*linedata(:,6);       totalbuses = max(max(linedata(:,1)),max(linedata(:,2)));    % total buses totalbranches = length(linedata(:,1));                      % no. of branches ybus = zeros(totalbuses,totalbuses);   for b=1:totalbranches     ybus((linedata(b,1)),(linedata(b,2)))=-y(b);     ybus((linedata(b,2)),(linedata(b,1))) =ybus((linedata(b,1)),(linedata(b,2)));      end  for c=1:totalbuses      for d=1:totalbranches          if linedata(d,1) == c || linedata(d,2) == c              ybus(c,c) = ybus(c,c) + y(d) + ss(d);          end      end  end disp('TABLE 9.3 PAGE # 338   BUS ADMITTANCE MATRIX FOR EXAMPLE 9.2') ybus z=zeros(totalbuses,4); busnumber=busdata(:,1); PG=busdata(:,2); QG=busdata(:,3); PL=busdata(:,4); QL=busdata(:,5); V=busdata(:,6); VV=V; ANG=busdata(:,7); type = busdata(:,8); P = (PG-PL)./100;    % per unit active power at buses Q = (QG-QL)./100;    % per unit reactive power at buses tol=1; iter=0; kk=0.0001; while tol > kk          for i = 2:totalbuses            YV = 0;                  for k = 1:totalbuses             if i~=k                 YV = YV + ybus(i,k)* V(k);  % multiplying admittance & voltage             YV;             end                      YV;         end         if busdata(i,8) == 3     %Calculating Qi for PV bus             %Q(i) = -imag(conj(V(i))*(YV + ybus(i,i)*V(i)));             Q(i) =  -imag(conj(V(i))*(YV + ybus(i,i)*V(i)));             busdata(i,3)=Q(i);         end                   % end         V(i) = (1/ybus(i,i))*((P(i)-j*Q(i))/conj(V(i)) - YV); % Compute Bus Voltages.                   % Calculating Corrected Voltage for PV bus        if busdata(i,8) == 3                    vc(i)=abs(VV(i))*(V(i)/abs(V(i)));            busdata(i,6)=vc(i);            V(i)=vc(i);        end                   % Calculating Accelerated Voltage for PQ bus        if busdata(i,8) == 2         VACC(i)= VV(i)+(V(i)-VV(i));         busdata(i,6)=VACC(i);         V(i)=VACC(i);        end         %V(i)=V;             end       iter = iter + 1;      % Increment iteration count.     tol = max(abs(abs(V) - abs(VV)));     % Calculate tolerance.    VV = V; end Q; iter; YV; V; %real(VACC') z(1:totalbuses,1)=busdata(:,1); z(1:totalbuses,2)=busdata(:,8); z(1:totalbuses,3)=abs(busdata(:,6)); z(1:totalbuses,4)=radtodeg(angle(V)); disp('           Bus No.   Bus Type     Voltage      Angle      '); z

hp prime porting code ...tried but got different results
Code:
 #cas gausss(M1,M3):= BEGIN  //M1=linedata, M2=ybus  M3= busdata ss:= .* col(M1,5);    y:= col(M1,3); y:=1 ./ y; totalbuses:= MAX(MAX(col(M1,1)),MAX(col(M1,2))); totalbranches:= length(M1); M2:=makemat(totalbuses,totalbuses); FOR b FROM 1 TO totalbranches DO     M2(M1(b,1),M1(b,2)):=-y(b);     M2((M1(b,2)),(M1(b,1))):=M2((M1(b,1)),(M1(b,2)));    END; FOR c FROM 1 TO totalbuses DO     FOR d FROM 1 TO totalbranches DO              IF M1(d,1)==c OR M1(d,2) == c THEN                 M2(c,c) := M2(c,c) + y(d)+ss(d);              END;     END; END; M2; M9:=makemat(totalbuses,4); busnumber:=col(M3,1); pg:=col(M3,2); qg:=col(M3,3); pl:=col(M3,4); ql:=col(M3,5); v:= col(M3,6); vv:=v; ANG:=  col(M3,7); types:= col(M3,8); p := (pg - pl) ./100; q := (qg - ql) ./100; toll:=1; iter:=0; kk:=0.0001 alfa:=1.6; WHILE toll > kk DO          FOR j FROM 2 TO totalbuses DO            YV := 0;                  FOR k FROM 1 TO totalbuses DO             IF j≠k THEN                 YV := YV + M2(j,k)* v(k);              END;         END; //Calculating Qi for PV bus         IF M3(j,8) == 2 THEN              q(j) :=  -IM((CONJ(v(j)))*(YV + M2(j,j)*v(j)));             M3(j,3):=q(j);         END;           v(j) := (1/M2(j,j))*((p(j)-*q(j))/CONJ(v(j)) - YV); //Compute Bus Voltages.          //Calculating Corrected Voltage for PV bus               IF M3(j,8) == 3 THEN                    vc(j):=(abs(vv(j)))*(v(j)/abs(v(j)));            M3(j,6):=vc(j);            v(j):=vc(j);        END; //Calculating Accelerated Voltage for PQ bus        IF M3(j,8) == 2 THEN           VACC(j):= vv(j)+alfa*(v(j)-vv(j));           M3(j,6):=VACC(j);           v(j):=VACC(j);        END;         END;      iter := iter + 1;   toll :=abs(abs(MAX(v))-abs(MAX(vv)));    vv := v; END; q; YV; v; END; #end

Attached File(s) Thumbnail(s)

 « Next Oldest | Next Newest »

 Messages In This Thread Help me port this Matlab code......pls....(ported)[Solved] - toshk - 11-25-2015 05:16 PM RE: Help me port this Matlab code......pls - Marcio - 11-25-2015, 06:46 PM RE: Help me port this Matlab code......pls - toshk - 11-25-2015, 07:06 PM RE: Help me port this Matlab code......pls - Marcio - 11-25-2015, 07:30 PM Power-Flow.....Gauss-Siedel Method - toshk - 11-29-2015, 07:18 AM

User(s) browsing this thread: 1 Guest(s)