While examining current sources, I figured that since the basis for the circuit solver is mesh analysis, the current source will have to be converted to a voltage source in series with a resistance. This means a voltage source and a resistance in series connected in parallel with a voltmeter. So essentially the voltage source has to be adjusted with respect to the voltage measured so as to produce the necessary current.
So this starts another task that I was planning to do just before I released a major version. That is of user defined objects. Suppose, a user wants to create an object which is a voltage source in series with an RLC network and make a dozen copies of it in the circuit, how would that be done? As of now, the entire solver works on the basis of cell positions on the spreadsheet. But now these cell positions won't exist in real but fictitious cell positions will have to be created.
Suppose such an object is to be defined. How would the details be internally? Suppose it is placed at position 7E in the spreadsheet. How would each element be defined as cell positions in the loop as it may be more than one loop if there is a parallel element. One possibility is to go ahead with just one cell position but replicate the loop by the number of parallel elements in the loop (essentially the number of branches in parallel). So every time the object is called, an identifier needs to be created to know which parallel element to consider.
So essentially a third dimension needs to be added. This makes sense as we are adding layers to the spreadsheet. But this might need a complete rework of the entire solver! Need to think over this.
Friday, April 26, 2013
Monday, April 22, 2013
Capacitor class
Been a long time since I updated the blog.
To add another element to the library - the capacitor. Anyway, the code is just an extension/combination of the voltmeter/voltage source classes. Here is the code (click on "view raw" below the code box to see it in a new window):
I'll release another version after completing the current source class which I hope will be in a few days.
For the current source class, my guess it would be best to exclude it from the branch list so that it doesn't appear in any of the loops. For every current source a simple nodal equation can then be added.
To add another element to the library - the capacitor. Anyway, the code is just an extension/combination of the voltmeter/voltage source classes. Here is the code (click on "view raw" below the code box to see it in a new window):
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| class Capacitor: | |
| """ Capacitor class. Contains functions to initiliaze | |
| the resistor according to name tag, unique cell position, | |
| update system matrix on each iteration. """ | |
| def __init__(self, cap_index, cap_pos, cap_tag): | |
| """ Constructor to initialize value. | |
| Also, takes in the identifiers - | |
| index (serial number), cell position and tag. """ | |
| self.type="Capacitor" | |
| self.cap_number=cap_index | |
| self.cap_pos=cap_pos | |
| self.cap_tag=cap_tag | |
| self.capacitor=10.0e-6 | |
| self.current=0.0 | |
| self.voltage=0.0 | |
| self.v_dbydt=0.0 | |
| self.cap_polrty=[-1, -1] | |
| def display(self): | |
| """ Displays info about the component.""" | |
| print "Capacitor is ", | |
| print self.cap_tag, | |
| print "= %f" %self.capacitor, | |
| print " located at ", | |
| print self.cap_pos, | |
| print " with positive polarity towards %s" %(csv_element(self.cap_polrty)) | |
| def ask_values(self, x_list, ckt_mat, sys_branch): | |
| """ Writes the values needed to the spreadsheet.""" | |
| cap_params=["Capacitor"] | |
| cap_params.append(self.cap_tag) | |
| cap_params.append(self.cap_pos) | |
| cap_params.append(self.capacitor) | |
| if self.cap_polrty==[-1, -1]: | |
| # Looking for a default value of polarity | |
| # in the neighbouring cells | |
| self.cap_elem=csv_tuple(self.cap_pos) | |
| if self.cap_elem[0]>0: | |
| if ckt_mat[self.cap_elem[0]-1][self.cap_elem[1]]: | |
| self.cap_polrty=[self.cap_elem[0]-1, self.cap_elem[1]] | |
| if self.cap_elem[1]>0: | |
| if ckt_mat[self.cap_elem[0]][self.cap_elem[1]-1]: | |
| self.cap_polrty=[self.cap_elem[0], self.cap_elem[1]-1] | |
| if self.cap_elem[0]<len(ckt_mat)-1: | |
| if ckt_mat[self.cap_elem[0]+1][self.cap_elem[1]]: | |
| self.cap_polrty=[self.cap_elem[0]+1, self.cap_elem[1]] | |
| if self.cap_elem[1]<len(ckt_mat)-1: | |
| if ckt_mat[self.cap_elem[0]][self.cap_elem[1]+1]: | |
| self.cap_polrty=[self.cap_elem[0], self.cap_elem[1]+1] | |
| else: | |
| for c1 in range(len(sys_branch)): | |
| if csv_tuple(self.cap_pos) in sys_branch[c1]: | |
| if not self.cap_polrty in sys_branch[c1]: | |
| print "!"*50 | |
| print "ERROR!!! Capacitor polarity should be in the same branch as the capacitor. Check source at %s" %self.cap_pos | |
| print "!"*50 | |
| cap_params.append("Positive polarity towards (cell) = %s" %csv_element(self.cap_polrty)) | |
| x_list.append(cap_params) | |
| def get_values(self, x_list, ckt_mat): | |
| """ Takes the parameter from the spreadsheet.""" | |
| self.capacitor=float(x_list[0]) | |
| cap_polrty=x_list[1].split("=")[1] | |
| # Convert the human readable form of cell | |
| # to [row, column] form | |
| while cap_polrty[0]==" ": | |
| cap_polrty=cap_polrty[1:] | |
| self.cap_polrty=csv_tuple(cap_polrty) | |
| if not ckt_mat[self.cap_polrty[0]][self.cap_polrty[1]]: | |
| print "Polarity incorrect or changed. Branch does not exist at %s" %csv_element(self.cap_polrty) | |
| def transfer_to_sys(self, sys_loops, mat_e, mat_a, mat_b, mat_u, source_list): | |
| """ The matrix B in E.dx/dt=Ax+Bu will be updated by the | |
| polarity of the capacitor.""" | |
| for c1 in range(len(sys_loops)): | |
| for c2 in range(len(sys_loops[c1][c1])): | |
| if csv_tuple(self.cap_pos) in sys_loops[c1][c1][c2]: | |
| # If the positive polarity appears before the capacitor position | |
| # it means as per KVL, we are moving from +ve to -ve | |
| # and so the capacitor voltage will be taken negative | |
| if sys_loops[c1][c1][c2].index(self.cap_polrty)<sys_loops[c1][c1][c2].index(csv_tuple(self.cap_pos)): | |
| if sys_loops[c1][c1][c2][-1]=="forward": | |
| mat_b.data[c1][source_list.index(self.cap_pos)]=-1.0 | |
| else: | |
| mat_b.data[c1][source_list.index(self.cap_pos)]=1.0 | |
| else: | |
| if sys_loops[c1][c1][c2][-1]=="forward": | |
| mat_b.data[c1][source_list.index(self.cap_pos)]=1.0 | |
| else: | |
| mat_b.data[c1][source_list.index(self.cap_pos)]=-1.0 | |
| return | |
| def transfer_to_branch(self, sys_branch, source_list): | |
| """ Transfers parameters to system branch if capacitor | |
| exists in the branch. """ | |
| if csv_tuple(self.cap_pos) in sys_branch: | |
| if sys_branch.index(self.cap_polrty)<sys_branch.index(csv_tuple(self.cap_pos)): | |
| sys_branch[-1][1][source_list.index(self.cap_pos)]=-1.0 | |
| else: | |
| sys_branch[-1][1][source_list.index(self.cap_pos)]=1.0 | |
| return | |
| def generate_val(self, source_lst, sys_loops, mat_u, t, dt): | |
| """ The capacitor voltage is updated in the matrix u in | |
| E.dx/dt=Ax+Bu .""" | |
| self.v_dbydt=self.current/self.capacitor | |
| self.voltage+=self.v_dbydt*dt | |
| mat_u.data[source_lst.index(self.cap_pos)][0]=self.voltage | |
| self.op_value=self.voltage | |
| def update_val(self, sys_loops, lbyr_ratio, mat_e, mat_a, mat_b, state_vec, mat_u): | |
| """ The capacitor current is calculated as a result of the KVL.""" | |
| self.current=0.0 | |
| for c1 in range(len(sys_loops)): | |
| for c2 in range(len(sys_loops[c1][c1])): | |
| if csv_tuple(self.cap_pos) in sys_loops[c1][c1][c2]: | |
| # If capacitor polarity is before the capacitor | |
| # position, it means current is positive. | |
| if sys_loops[c1][c1][c2].index(self.cap_polrty)<sys_loops[c1][c1][c2].index(csv_tuple(self.cap_pos)): | |
| # Then check is the loop is aiding or opposing | |
| # the main loop. | |
| if sys_loops[c1][c1][c2][-1]=="forward": | |
| self.current+=state_vec.data[c1][0] | |
| else: | |
| self.current-=state_vec.data[c1][0] | |
| else: | |
| if sys_loops[c1][c1][c2][-1]=="forward": | |
| self.current-=state_vec.data[c1][0] | |
| else: | |
| self.current+=state_vec.data[c1][0] | |
| return |
I'll release another version after completing the current source class which I hope will be in a few days.
For the current source class, my guess it would be best to exclude it from the branch list so that it doesn't appear in any of the loops. For every current source a simple nodal equation can then be added.
Friday, April 5, 2013
Version 0.1.4 released
Released another version:
http://sourceforge.net/projects/pythonpowerelec/files/
Next step:
Complete the passive library - include capacitor and current source. Current source is extremely important for us power electronics guys and that might be a bit of a challenge as the loops will be refined every time a current source is found.
http://sourceforge.net/projects/pythonpowerelec/files/
Next step:
Complete the passive library - include capacitor and current source. Current source is extremely important for us power electronics guys and that might be a bit of a challenge as the loops will be refined every time a current source is found.
Simulation speed
Before the topic of simulation speed, a minor check to be performed on components that have polarity like voltage sources, ammeters and voltmeters.
Until now, the only check on the polarity of these elements is that it is within the physical circuit map. But now that branches have been defined, another check can be enforced about whether the polarity element is in the same branch as the component. This is the code (click on "view raw" below the box to see the code in a new window):
For this to happen, the main program "circuit_solver.py" has been rearranged a bit. No changes to the code as such. Just that the network interpreter is called immediately after the circuit layout file is interpreted and checked for basic errors like componenents not found in library or duplicate components. The idea is that polarity errors should be picked up at the design stage before the simulation starts.
Now the other change of simulation speed. My concept has been similar to an interrupt driver program. Instead of an interrupt vector, generate an "event vector". This event vector will be "yes" corresponding to a branch is an element has changed in the branch. Otherwise, it will remain "no". If there is a single "yes", this would mean the all the system matrices A, B and E would be calculated completely. Currently, this might be necessary. At a later stage, I could think of ways to recalculate only the parts that are affected.
So, the branch_events vector is reset to "no" at the end of the code block. This vector is initialized to "yes" at the beginning to ensure that at least one computation takes place. With passive circuits, there won't be any event generation and therefore, there will be no recalculation.
This has speeded up the execution quite a bit but still not as much as I would have liked. Need to figure which of the code blocks can be done away with. A fairly simple circuit like "testckt7.csv" in the example list which results in a 4x4 system takes around 2 minutes to solve for a run time of 0.5s and a time step of 10 microseconds. Which is quite slow. At later stages, this might be a problem when simulating an inverter with a time step of 1 microsecond. If the inverter has a switching frequency of 5 kHz, an event would be generated around every 30 microseconds. No idea how it is going to be.
Until now, the only check on the polarity of these elements is that it is within the physical circuit map. But now that branches have been defined, another check can be enforced about whether the polarity element is in the same branch as the component. This is the code (click on "view raw" below the box to see the code in a new window):
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| def ask_values(self, x_list, ckt_mat, sys_branch): | |
| """ Writes the values needed to the spreadsheet.""" | |
| volt_params=["VoltageSource"] | |
| volt_params.append(self.volt_tag) | |
| volt_params.append(self.volt_pos) | |
| volt_params.append("Peak (Volts) = %f" %self.v_peak) | |
| volt_params.append("Frequency (Hertz) = %f" %self.v_freq) | |
| volt_params.append("Phase (degrees) = %f" %self.v_phase) | |
| if self.v_polrty==[-1, -1]: | |
| # Looking for a default value of polarity | |
| # in the neighbouring cells | |
| self.volt_elem=csv_tuple(self.volt_pos) | |
| if self.volt_elem[0]>0: | |
| if ckt_mat[self.volt_elem[0]-1][self.volt_elem[1]]: | |
| self.v_polrty=[self.volt_elem[0]-1, self.volt_elem[1]] | |
| if self.volt_elem[1]>0: | |
| if ckt_mat[self.volt_elem[0]][self.volt_elem[1]-1]: | |
| self.v_polrty=[self.volt_elem[0], self.volt_elem[1]-1] | |
| if self.volt_elem[0]<len(ckt_mat)-1: | |
| if ckt_mat[self.volt_elem[0]+1][self.volt_elem[1]]: | |
| self.v_polrty=[self.volt_elem[0]+1, self.volt_elem[1]] | |
| if self.volt_elem[1]<len(ckt_mat)-1: | |
| if ckt_mat[self.volt_elem[0]][self.volt_elem[1]+1]: | |
| self.v_polrty=[self.volt_elem[0], self.volt_elem[1]+1] | |
| else: | |
| for c1 in range(len(sys_branch)): | |
| if csv_tuple(self.volt_pos) in sys_branch[c1]: | |
| if not self.v_polrty in sys_branch[c1]: | |
| print "!"*50 | |
| print "ERROR!!! Voltage source polarity should be in the same branch as the voltage source. Check source at %s" %self.volt_pos | |
| print "!"*50 | |
| volt_params.append("Positive polarity towards (cell) = %s" %csv_element(self.v_polrty)) | |
| x_list.append(volt_params) |
For this to happen, the main program "circuit_solver.py" has been rearranged a bit. No changes to the code as such. Just that the network interpreter is called immediately after the circuit layout file is interpreted and checked for basic errors like componenents not found in library or duplicate components. The idea is that polarity errors should be picked up at the design stage before the simulation starts.
Now the other change of simulation speed. My concept has been similar to an interrupt driver program. Instead of an interrupt vector, generate an "event vector". This event vector will be "yes" corresponding to a branch is an element has changed in the branch. Otherwise, it will remain "no". If there is a single "yes", this would mean the all the system matrices A, B and E would be calculated completely. Currently, this might be necessary. At a later stage, I could think of ways to recalculate only the parts that are affected.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # Check is an event has been generated. If yes, | |
| # recalculate the system matrices. If not, just | |
| # solve the ODE | |
| if "yes" in branch_events: | |
| # Make a copy of the system_loops | |
| # as row operations may be performed on them. | |
| system_loops_copy=[] | |
| for c1 in range(system_size): | |
| row_vector=[] | |
| for c2 in range(system_size): | |
| element_vector=[] | |
| for c3 in range(len(system_loops[c1][c2])): | |
| inter_vector=[] | |
| for c4 in range(len(system_loops[c1][c2][c3])): | |
| inter_vector.append(system_loops[c1][c2][c3][c4]) | |
| element_vector.append(inter_vector) | |
| row_vector.append(element_vector) | |
| system_loops_copy.append(row_vector) | |
| # Initialize the entries in branch_params | |
| for c1 in range(len(branch_params)): | |
| branch_params[c1][-1][0][0]=0.0 | |
| branch_params[c1][-1][0][1]=0.0 | |
| for c2 in range(len(source_list)): | |
| branch_params[c1][-1][1][c2]=0.0 | |
| # Initialize A,E and B matrices to zero. | |
| sys_mat_a.zeros(system_size, system_size) | |
| sys_mat_e.zeros(system_size, system_size) | |
| sys_mat_b.zeros(system_size, source_size) | |
| # To update the branch_params list. | |
| # Take every element in a branch and check | |
| # if that is a valid object identifier | |
| # If not, i.e if it is a wire, an exception will | |
| # be raised in which case nothing will happen. | |
| # If no exception is raised, do the transfer to the branch | |
| for c1 in range(len(branch_params)): | |
| for c2 in range(len(branch_params[c1][:-1])): | |
| try: | |
| comp_pos=csv_element(branch_params[c1][c2]) | |
| component_objects[comp_pos] | |
| except: | |
| pass | |
| else: | |
| component_objects[comp_pos].transfer_to_branch(branch_params[c1], source_list) | |
| # Now transfer values from the branch to the matrices A,B and E. | |
| # Look where each branch appears and transfer the info. | |
| for c1 in range(len(system_loops_copy)): | |
| for c2 in range(len(system_loops_copy)): | |
| for c3 in range(len(system_loops_copy[c1][c2])): | |
| for c4 in range(len(branch_params)): | |
| if system_loops_copy[c1][c2][c3][:-1]==branch_params[c4][:-1]: | |
| if c1==c2: | |
| sys_mat_a.data[c1][c2]+=branch_params[c4][-1][0][0] | |
| sys_mat_e.data[c1][c2]+=branch_params[c4][-1][0][1] | |
| if system_loops_copy[c1][c2][c3][-1]=="forward": | |
| for c5 in range(len(source_list)): | |
| sys_mat_b.data[c1][c5]+=branch_params[c4][-1][1][c5] | |
| else: | |
| for c5 in range(len(source_list)): | |
| sys_mat_b.data[c1][c5]-=branch_params[c4][-1][1][c5] | |
| else: | |
| if system_loops_copy[c1][c2][c3][-1]=="forward": | |
| sys_mat_a.data[c1][c2]+=branch_params[c4][-1][0][0] | |
| sys_mat_e.data[c1][c2]+=branch_params[c4][-1][0][1] | |
| else: | |
| sys_mat_a.data[c1][c2]-=branch_params[c4][-1][0][0] | |
| sys_mat_e.data[c1][c2]-=branch_params[c4][-1][0][1] | |
| # Check if the system is stiff and recalculate if it is. | |
| remove_stiffness(sys_mat_e, sys_mat_a, sys_mat_b, dt, system_loops_copy, branch_params) | |
| mat_ode_reduce(sys_mat_e, sys_mat_a, sys_mat_b) | |
| # Remove the branch events | |
| for c1 in range(len(branch_events)): | |
| branch_events[c1]="no" | |
| # Update the input values in u matrix | |
| for c1 in range(len(source_list)): | |
| component_objects[source_list[c1]].generate_val(source_list, system_loops_copy, sys_mat_u, t, dt) | |
| # The ODE solver. | |
| # Will return the x(k+1) value called | |
| # as next_state_vec from x(k) value | |
| # called curr_state_vec | |
| mat_ode(sys_mat_e, sys_mat_a, sys_mat_b, [curr_state_vec, next_state_vec], sys_mat_u, dt, ode_var, system_loops_copy, node_list) |
So, the branch_events vector is reset to "no" at the end of the code block. This vector is initialized to "yes" at the beginning to ensure that at least one computation takes place. With passive circuits, there won't be any event generation and therefore, there will be no recalculation.
This has speeded up the execution quite a bit but still not as much as I would have liked. Need to figure which of the code blocks can be done away with. A fairly simple circuit like "testckt7.csv" in the example list which results in a 4x4 system takes around 2 minutes to solve for a run time of 0.5s and a time step of 10 microseconds. Which is quite slow. At later stages, this might be a problem when simulating an inverter with a time step of 1 microsecond. If the inverter has a switching frequency of 5 kHz, an event would be generated around every 30 microseconds. No idea how it is going to be.
Tuesday, April 2, 2013
Version 0.1.3 released
Been incredibly busy with my official research project. Will continue to be busy for another week until my PCBs are sent out for fabrication. Thought I would release another version just so that I don't loose track of what I have done.
http://sourceforge.net/projects/pythonpowerelec/
Created a voltmeter model that you can use to measure voltage. The concept is it is a large resistance is parallel with the nodes across which voltage is to be measured. The parameter taken as input is the rated voltage of the system. This is just an approximate value so as to choose the voltmeter resistance large enough that a negligible current of 1 microampere flows through the voltmeter.
http://sourceforge.net/projects/pythonpowerelec/
Created a voltmeter model that you can use to measure voltage. The concept is it is a large resistance is parallel with the nodes across which voltage is to be measured. The parameter taken as input is the rated voltage of the system. This is just an approximate value so as to choose the voltmeter resistance large enough that a negligible current of 1 microampere flows through the voltmeter.
Subscribe to:
Posts (Atom)