+- 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 ');
zhp 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;
#endRE: 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