Nodal analysis

There might have been a mistake in the previous current initialization of branch currents. The branches are current sources if:
1. They have inductors and are not stiff.
2. They have voltages and have zero impedance.

If the branches are stiff or are purely resistive with or without a voltage, their currents will need to be calculated by nodal analysis. The code to set the branch currents is:

	# Initializing the branch currents for the next set of nodal analysis
	for c1 in range(len(branch_params)):
	branch_currents[c1]=0.0
	for c1 in range(len(branch_params)):
	# Check if the branch is stiff
	if stiff_ratio[c1]=="no":
	# Check if the branch has inductance
	if branch_params[c1][-1][0][1]:
	branch_currents[c1]=branch_params[c1][-1][2]
	# Check if it is a zero impedance branch with a voltage
	elif ((branch_params[c1][-1][0][0]==0.0) and (1.0 in branch_params[c1][-1][1] or -1.0 in branch_params[c1][-1][1])):
	branch_currents[c1]=branch_params[c1][-1][2]
	# This second run of the function is to determine the new
	# branch currents because of any possible change in the state
	# of nonlinear devices.
	current_continuity(node_list, branch_params, stiff_ratio, node_voltage, branch_currents, admittance_matrix, source_vector, sys_mat_u)

view raw nodal_current_init.py hosted with ❤ by GitHub

The actual nodal analysis is:

	# Start iterating through the nodes
	for c1 in range(len(list_of_nodes)):
	# Iterate through all the branches.
	for c2 in range(len(branch_info)):
	# Look for the node as the starting or ending node in each branche
	if ((list_of_nodes[c1]==branch_info[c2][0]) or (list_of_nodes[c1]==branch_info[c2][-2])):
	if (list_of_nodes[c1]==branch_info[c2][0]):
	start_node_pos=list_of_nodes.index(branch_info[c2][0])
	# Look for the start node in the short node list
	# If it is a short node, the first short node
	# in that list to which the current short node
	# belongs will be the start node pos.
	# this is because all short nodes are the same.
	for c3 in range(len(shortnode_list)):
	if start_node_pos in shortnode_list[c3]:
	start_node_pos=shortnode_list[c3][0]
	# Mark the position of the other node
	end_node_pos=list_of_nodes.index(branch_info[c2][-2])
	# Again check if it is a short node.
	for c3 in range(len(shortnode_list)):
	if end_node_pos in shortnode_list[c3]:
	end_node_pos=shortnode_list[c3][0]
	# The default direction of current is taken to be
	# away from the starting node - start to end.
	branch_src_dir=1.0
	else:
	start_node_pos=list_of_nodes.index(branch_info[c2][-2])
	for c3 in range(len(shortnode_list)):
	if start_node_pos in shortnode_list[c3]:
	start_node_pos=shortnode_list[c3][0]
	end_node_pos=list_of_nodes.index(branch_info[c2][0])
	for c3 in range(len(shortnode_list)):
	if end_node_pos in shortnode_list[c3]:
	end_node_pos=shortnode_list[c3][0]
	branch_src_dir=-1.0
	# If resistance of the branch is non-zero
	# if branch_info[c2][-1][0][0]:
	# # If the branch has no inductance or is a stiff branch
	# # in which case inductance is negligible compared
	# # to the resistance and the branch is stiff.
	# if (branch_info[c2][-1][0][1]==0.0):
	# mho_matrix[start_node_pos][start_node_pos]+=1.0/branch_info[c2][-1][0][0]
	# mho_matrix[start_node_pos][end_node_pos]-=1.0/branch_info[c2][-1][0][0]
	#
	# # Similarly, any voltage source that exists
	# # will be treated as positive if it forces a
	# # a current away the starting node.
	# # The way branch_info contains these sources is
	# # compatible to the definition
	# # src_vector is on the RHS of the equation,
	# # so thus the -ve sign
	# for c3 in range(len(branch_info[c2][-1][1])):
	# src_vector[start_node_pos]-=branch_src_dir*branch_info[c2][-1][1][c3]*sys_inputs.data[c3][0]/branch_info[c2][-1][0][0]
	# # Check if inductance is non-zero
	# if branch_info[c2][-1][0][1]:
	# # Check if branch is not stiff
	# if branch_stiff[c2]=="no":
	# # Add the current to the src_vector as a current source
	# src_vector[start_node_pos]-=branch_src_dir*br_currents[c2]
	# # Check if is a zero impedance branch
	# if ((branch_info[c2][-1][0][0]==0.0) and (branch_info[c2][-1][0][1]==0.0)):
	# if ((1.0 in branch_info[c2][-1][1]) or (-1.0 in branch_info[c2][-1][1])):
	# src_vector[start_node_pos]-=branch_src_dir*br_currents[c2]
	# Check if branch is not stiff
	if branch_stiff[c2]=="no":
	# Check if inductance is non-zero
	if branch_info[c2][-1][0][1]:
	src_vector[start_node_pos]-=branch_src_dir*br_currents[c2]
	# Check if is a zero impedance branch with a voltage
	elif ((branch_info[c2][-1][0][0]==0.0) and (1.0 in branch_info[c2][-1][1] or -1.0 in branch_info[c2][-1][1])):
	src_vector[start_node_pos]-=branch_src_dir*br_currents[c2]
	# Check if it is a resistive branch
	if ((branch_info[c2][-1][0][0]) and (branch_info[c2][-1][0][1]==0.0)):
	mho_matrix[start_node_pos][start_node_pos]+=1.0/branch_info[c2][-1][0][0]
	mho_matrix[start_node_pos][end_node_pos]-=1.0/branch_info[c2][-1][0][0]
	# Similarly, any voltage source that exists
	# will be treated as positive if it forces a
	# a current away the starting node.
	# The way branch_info contains these sources is
	# compatible to the definition
	# src_vector is on the RHS of the equation,
	# so thus the -ve sign
	for c3 in range(len(branch_info[c2][-1][1])):
	src_vector[start_node_pos]+=branch_src_dir*branch_info[c2][-1][1][c3]*sys_inputs.data[c3][0]/branch_info[c2][-1][0][0]

view raw nodal_analysis.py hosted with ❤ by GitHub

Other than that, the code is quite the same as the previous case. This code was tested with a three phase diode bridge rectifier and was found to be working OK. The previously released code was breaking.