Help me port this Matlab code......pls....(ported)[Solved] - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html) +--- Forum: HP Prime (/forum-5.html) +--- Thread: Help me port this Matlab code......pls....(ported)[Solved] (/thread-5212.html) Help me port this Matlab code......pls....(ported)[Solved] - toshk - 11-25-2015 05:16 PM 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``` RE: Help me port this Matlab code......pls - Marcio - 11-25-2015 06:46 PM 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? RE: Help me port this Matlab code......pls - toshk - 11-25-2015 07:06 PM (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 [attachment=2864] RE: Help me port this Matlab code......pls - Marcio - 11-25-2015 07:30 PM (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. Power-Flow.....Gauss-Siedel Method - toshk - 11-29-2015 07:18 AM 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```