A few month ago, mystran published on KVR a small SPICE simulator for real-time processing. I liked the idea, the drawback being that the code is generic and not tailored like a static version of the optimizer. So I wondered if it was doable. But for this, I have to start from the basics and build from there. So let’s go.

Why Python?

I’ve decided to do the prototype in Python. The reason is simple, it’s easy to create something very fast in Python, write all the basic tests there and figure out what functionalities are required.

First, the objective is to have a statically generated model in the long-term, so I need to differentiate between static voltage pins, input voltage pins and output voltage pins. The latter ones could be any pins for which the voltage will be computed by the modeller.

Then we need components, like resistors, diodes, transistors… Capacitors and coils are also required, but let’s use the model I presented in a previous blog post. This will simplify writing the equations.

So what is the basic equation in an electronic modeller? The easiest approach is nodal analysis, something small and easy to understand, like Kirchhoff’s current law . It is simple enough and (almost) all the models I use describe current from a pin as a function of voltages.

There is one thing in the prototype that is not perfect and that I haven’t figured out yet. A circuit is in a steady state when we start feeding it an audio signal. To compute it, we need to make capacitor like open circuits (easy) and coil like short circuits (not easy). The issue with short-circuits also happens when you have a variable resistor than you turn off entirely. In a dynamic model, we can easily collapse pins together when required, but in a static model, when you want to optimize the shape of the matrices, this is less than ideal. So to avoid this, I use a very small resistor. Not perfect, but it seems to work. For now.

Description of the different methods

Let’s start with the basic Modeler class (I wrote the prototype with American English and the C++ version in British English, seems like I can’t decide which one I should use…). The constructor takes a number of dynamic pins (the ones we will compute the voltage for), static pins (fixed input voltage) and input pins (variable input voltage). I will keep a list of the components of the model, then a structure for accessing the pins (useful for the dynamic pins so that we can get the sum of the currents) and the same for the state voltages. The distinction between dynamic, static and input voltages will be done through the character ‘D’, ‘S’ and ‘I’ (in Python, in C++, I’ll use an enum). We actually don’t need to store anything more than the dynamic pins, but for the sake of this prototype, let’s store all of them.

class Modeler(object):
    """
    Modeling class
    """
    def __init__(self, nb_dynamic_pins, nb_static_pins, nb_inputs = 0):
        self.components = []
        self.dynamic_pins = [[] for i in range(nb_dynamic_pins)]
        self.static_pins = [[] for i in range(nb_static_pins)]
        self.input_pins = [[] for i in range(nb_inputs)]
        self.pins = {
                'D': self.dynamic_pins,
                'S': self.static_pins,
                'I': self.input_pins,
                }
 
        self.dynamic_state = np.zeros(nb_dynamic_pins, dtype=np.float64)
        self.static_state = np.zeros(nb_static_pins, dtype=np.float64)
        self.input_state = np.zeros(nb_inputs, dtype=np.float64)
        self.state = {
                'D': self.dynamic_state,
                'S': self.static_state,
                'I': self.input_state,
                }
        self.initialized = False

To add a new component, we need to keep the component and store the pins inside the component (because the components needs to track the pins for the current computation) and for each pin, we store the component attached to it and the index of that pin for the component. The latter functionality will allow to get the proper sign of the current for the pin when we compute our equations.

    def add_component(self, component, pins):
        """
        add a new component
        :param component: component to add
        :param pins: list of tuples indicating how the component is connected
        """
        self.components.append(component)
        component.pins = pins
        for (i, pin) in enumerate(pins):
            t, pos = pin
            self.pins[t][pos].append((component, i))

As I’ve said, we need to compute a steady state. When we do so, we need to start by initializing components (coils and capacitors), solve for a steady state and then update the same components to update their internal state.

    def setup(self):
        """
        Initializes the internal state
        :param steady_state: if set to True (default), computes a steady state
        """
        for component in self.components:
            component.update_steady_state(self.state, self.dt)
 
        self.solve(True)
 
        for component in self.components:
            component.update_steady_state(self.state, self.dt)
 
        self.initialized = True

To update the state based on an input, we do the following, reusing our solve() method:

    def __call__(self, input):
        """
        Works out the value for the new input vector
        :param input: vector of input values
        """
        if not self.initialized:
            self.setup()
 
        self.input_state[:] = input
 
        self.solve(False)
 
        for component in self.components:
            component.update_state(self.state)
 
        return self.dynamic_state

Now we can write the solver part. We iterate several time, and for each component, we tell them to precompute their state (for costly components like diodes, transistors or valves), and then for each pin, we write an equation (remember that each pin is a state we have to solve, so with as many equations as we have pins, we should be able to get to the next state) and get the Jacobian for that equation. If the Kirchhoff equations are already satisfied (i.e. close to 0), we stop. If the delta we compute is small enough, we also stop.

    def solve(self, steady_state):
        """
        Actually solve the equation system
        :param steady_state: if set to True (default), computes for a steady state
        """
        iteration = 0
        while iteration < MAX_ITER and not self.iterate(steady_state):
            iteration = iteration + 1        
 
    def iterate(self, steady_state):
        """
        Do one iteration
        :param steady_state: if set to True (default), computes for a steady state
        """
        for component in self.components:
            component.precompute(self.state, steady_state)
        eqs = []
        jacobian = []
        for i, pin in enumerate(self.dynamic_pins):
            eq, jac = self.compute_current(pin, steady_state)
            eqs.append(eq)
            jacobian.append(jac)
 
        eqs = np.array(eqs)
        jacobian = np.array(jacobian)
 
        if np.all(np.abs(eqs) < EPS):
            return True
 
        delta = np.linalg.solve(jacobian, eqs)
        if np.all(np.abs(delta) < EPS):
            return True
 
        self.dynamic_state -= delta
 
        return False

Now the last missing part will be the actual equation and Jacobian line building from the components:

    def compute_current(self, pin, steady_state):
        """
        Compute Kirschhoff law for the non static pin
        Compute also the jacobian for all the connected pins
        :param pin: tuple indicating which pin we compute the current for
        :param steady_state: if set to True (default), computes for a steady state
        """
        eq = sum([component.get_current(i, self.state, steady_state) for (component, i) in pin])
        jac = [0] * len(self.dynamic_state)
        for (component, j) in pin:
            for (i, component_pin) in enumerate(component.pins):
                if component_pin[0] == "D":
                    jac[component_pin[1]] += component.get_gradient(j, i, self.state, steady_state)
        return eq, jac
 
def retrieve_voltage(state, pin):
    """
    Helper function to get the voltage for a given pin
    """
    return state[pin[0]][pin[1]]

Thanks to the way the modeller is built, we can pass the entire state and keep ‘D’, ‘S’ and ‘I’ to check for which voltage we need to compute the Jacobian.

Conclusion

As this blog post is already long, I’ll pause here for now before tackling the different components and what we need to change to the modeller for some ideal components like opamp.

The code is available on GitHub.