Demo entry 6661138

jacobi

   

Submitted by wf on Nov 20, 2017 at 11:18
Language: Matlab. Code size: 1.4 kB.

%% jacobi: function description
function [H,L,M,N] = jacobi(PQ_PV_mat,ybus,P,Q,nodes_num,PQ_num)
	H = zeros(nodes_num-1,nodes_num-1);
	L = zeros(PQ_num,PQ_num);
	M = zeros(PQ_num,nodes_num-1);
	N = zeros(nodes_num-1,PQ_num);
	for m = 1:nodes_num-1
		for n = 1:nodes_num-1
			Um = PQ_PV_mat(m,4);
			Un = PQ_PV_mat(n,4);
			Pm = P(m);
			Qm = Q(m);
			theta_mn = deg2rad(PQ_PV_mat(m,5)-PQ_PV_mat(n,5));
			Gmn = real(ybus(m,n));
			Bmn = imag(ybus(m,n));
			if m == n
				H(m,m) = Qm+Bmn*Um^2;
				%diagonal terms
				if n <= PQ_num
					N(m,m) = -Pm-Gmn*Um^2;
				end
			else
				H(m,n) = -Um*Un*(Gmn*sin(theta_mn)-Bmn*cos(theta_mn));
				%off diagronal terms
				if n <= PQ_num
					N(m,n) = -Um*Un*(Gmn*cos(theta_mn)+Bmn*sin(theta_mn));
				end
			end
		end
	end
	for m = 1:PQ_num
		for n = 1:nodes_num-1
			Um = PQ_PV_mat(m,4);
			Un = PQ_PV_mat(n,4); 
			Pm =P(m);
			Qm = Q(m);
			theta_mn = deg2rad(PQ_PV_mat(m,5)-PQ_PV_mat(n,5));
			Gmn = real(ybus(m,n));
			Bmn = imag(ybus(m,n));
			if m == n
				%m==n
				M(m,m) = -Pm+Gmn*Um^2;
				%diagonal terms
				if n <= PQ_num
					L(m,m) = -Qm+Bmn*Um^2;
				end
			else
				%m!=n
				M(m,n) = Um*Un*(Gmn*cos(theta_mn)+Bmn*sin(theta_mn));
				%off diagronal terms
				if n <= PQ_num
					L(m,n) = - Um*Un*(Gmn*sin(theta_mn)-Bmn*cos(theta_mn));
				end
			end
		end
	end

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).