Help me port this Matlab code......pls....(ported)[Solved]

[toshk]

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

[Marcio]

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?

[toshk]

(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
   

[Marcio]

(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.

[toshk]

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