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: 193 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)

11-25-2015, 06:46 PM
Post: #2
 Marcio Senior Member Posts: 438 Joined: Feb 2015
RE: Help me port this Matlab code......pls
I might be able to help, but first it'd be interesting see what comes up if you debug them simultaneously, line by line.

Does the Prime code have any syntax error?
11-25-2015, 07:06 PM (This post was last modified: 11-25-2015 07:12 PM by toshk.)
Post: #3
 toshk Member Posts: 193 Joined: Feb 2015
RE: Help me port this Matlab code......pls
(11-25-2015 06:46 PM)Marcio Wrote:  I might be able to help, but first it'd be interesting see what comes up if you debug them simultaneously, line by line.

Does the Prime code have any syntax error?

no major syntax error....
but the while loop fails to iterate.
suspect: % Calculate tolerance. in Matlab translation.
also some variants in YV values.
input data

11-25-2015, 07:30 PM
Post: #4
 Marcio Senior Member Posts: 438 Joined: Feb 2015
RE: Help me port this Matlab code......pls
(11-25-2015 07:06 PM)toshk Wrote:  no major syntax error....
but the while loop fails to iterate.

Eliminate all syntax errors and see if the loop finishes.
11-29-2015, 07:18 AM (This post was last modified: 11-30-2015 04:00 AM by toshk.)
Post: #5
 toshk Member Posts: 193 Joined: Feb 2015
Power-Flow.....Gauss-Siedel Method
Code:
 #cas gausss(M1,M3):= BEGIN //M1 for linedata: column1=bus From, column2=bus To, column3=Zpu, column4=Shunt Admittance [Q value] but can enter zero if not know , column5=Shunt Admittance [Bpu/2] //M3 for Busdata: column1=buses, column2=Pg, column3=Qg, column4=PL, column5=QL, column6=Vpu, column7= Vpu Angle, column8=Bus Type //Bus Type: 1.Slack Bus    2.PQ Bus    3.PV Bus ss:=(col(M1,5)) .* ; // capacitor    y:=(col(M1,3)); y:=approx(1 ./ y); totalbuses:= MAX(MAX(col(M1,1)),MAX(col(M1,2))); totalbranches:= length(col(M1,1)); 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;//ybus matrix busnumber:=col(M3,1); pg:=col(M3,2); qg:=col(M3,3); pl:=col(M3,4); ql:=col(M3,5); L2:= col(M3,6); vv:=L2; ANG:= col(M3,7); types:=col(M3,8); p:= approx((pg - pl) ./100); //Sbase of 100MVA L0:=approx((qg - ql) ./100); toll:=1; iter:=0; kk:=0.01;//tolerance 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)* L2(k);              END;         YV;         END; //Calculating Qi for PV bus         IF M3(j,8)== 3 THEN             L0(j) :=  -im(conj(L2(j))*(YV + M2(j,j)*L2(j)));             M3(j,3):=L0(j);         END;          v1:=L2(j);         m2:=M2(j,j);         pu:=p(j);         qu:=L0(j);         L2(j):= (1/m2)*((pu-*qu)/conj(v1) - YV); //Compute Bus Voltages.          //Calculating Corrected Voltage for PV bus               IF M3(j,8)== 3 THEN            L4(j):=abs(vv(j))*(L2(j)/abs(L2(j)));            M3(j,6):=L4(j);            L2(j):=L4(j);        END; //Calculating Accelerated Voltage for PQ bus        IF M3(j,8)== 2 THEN           L5(j):= vv(j)+alfa*(L2(j)-vv(j));           M3(j,6):=L5(j);           L2(j):=L5(j);        END;         END;      iter := iter + 1; L1:=vv;   toll:= max(abs(abs(L2)- abs(L1)));    vv:= L2; END; msgbox("check M3 col 6 for |V| and Angle; type M2 for Ybus ") END; #end
 « Next Oldest | Next Newest »

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