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