#Solving circuits with sympy

This documents explains how to solve electrical circuits using the **sympy** symbolic math module.

# Imports

First we need to import the sympy module


```python
# Import the sympy module
import sympy
```

## Example DC circuit

The following circuit includes one voltage source, one current source and two resistors.

The objective is to obtain the output voltage **Vo** as function of the components.



The circuit will be solved using the **nodal** method.

First we need to locate the circuit **nodes**, assign one as **ground** and assign a number for the rest of them. As the circuit has three nodes and one is ground, we have two nodes left: 1 and 2.

We will first generate a set of sympy symbols. There will be:

* One symbol for each component: Vs, R1, R2, Is

* One current symbol for each power sypply: iVs

* One symbol for each measurement we want to obtain: Vo

* One symbol for each node voltage that is not ground: V1, V2


```python
# Create the circuit symbols
Vs,iVs,R1,R2,Is,Vo, V1, V2 = sympy.symbols('Vs,iVs,R1,R2,Is,Vo,V1,V2')
```

Then we can define the current equations on each node except ground.

The current equations add all the currents in the node from each component.

All equations we add, are supposed to have a result of **zero**.


```python
# Create an empty list of equations
equations = []

# Nodal equations
equations.append(iVs-(V1-V2)/R1)       # Node 1
equations.append(Is-(V2-V1)/R1-V2/R2)  # Node 2
```

Then we add one equation for each voltage source that associates its voltage with the node voltages.

If we want to use two sided equations, we can use the sympy **Eq** function that equates the two sides.


```python
# Voltage source equations
equations.append(sympy.Eq(Vs,V1))
```

Finally we add one equation for each measurement we want to obtain


```python
# Measurement equations
equations.append(sympy.Eq(Vo,V2))
```

Now we can define the unknows for the circuit.

The number of unknowns shall be equal to the number of equations. 

The list includes:

* The node voltages: V1, V2

* The current on voltage sources: iVs

* The measurement values: Vo


```python
unknowns = [V1,V2,iVs,Vo]
```

We can see the equations and unknows before solving the circuit.

To ease reusing the code, we will define a **showCircuit** function that shows equations and unknowns


```python
# Define the function
def showCircuit():
    print('Equations')
    for eq in equations:
        print('    ',eq)
    print()
    print('Unknowns:',unknowns)
    print()
    
# Use the function
showCircuit()
```

    Equations
         iVs - (V1 - V2)/R1
         Is - V2/R2 - (-V1 + V2)/R1
         Eq(Vs, V1)
         Eq(Vo, V2)
    
    Unknowns: [V1, V2, iVs, Vo]
    


Now, we can solve the circuit.

The sympy **solve** function gets a list of equations and unknowns and return a **dictionary** with solved unknowns

The following code solves the circuit and list the solutions


```python
# Solve the circuit
solution = sympy.solve(equations,unknowns)

# List the solutions
print('Solutions')
for sol in solution:
    print('  ',sol,'=',solution[sol])
```

    Solutions
       V1 = Vs
       V2 = R2*(Is*R1 + Vs)/(R1 + R2)
       iVs = (-Is*R2 + Vs)/(R1 + R2)
       Vo = R2*(Is*R1 + Vs)/(R1 + R2)


Note that, in this case, the equation that includes **iVs** is only needed to obtain this unknown, so we can eliminate its equation and the **iVs** unknown if we don't need the **iVs** solution.

Note also that we can easily identify **V1** as **Vs**, so we can also eliminate the **Vs** equation if we use **Vs** as the voltage on node **V1**

Finally note that we can also identify **V2** as **Vo** so we can also eliminate the **Vo** equation.

The following code solves the circuit using only one equations


```python
solution = sympy.solve(Is-(Vo-Vs)/R1-Vo/R2,Vo)
print('Vo =',solution[0])
```

    Vo = R2*(Is*R1 + Vs)/(R1 + R2)


## Solve using the loop current method

Instead of using the **nodal** method, we could also use the complementary **loop current** method



In this method we assign a current to each loop in the circuit

We will first generate a set of sympy symbols. There will be:

* One symbol for each component: Vs, R1, R2, Is

* One voltage symbol for each current sypply: Vo

* One symbol for each measurement we want to obtain: Vo

* One symbol for each loop current: I1, I2

Note that, in this circuit, **Vo** appears two times, as the voltage on the **Is** source and as the value to measure. Logically we only define one **Vo** symbol.


```python
# Create the circuit symbols
Vs,R1,R2,Is,Vo,I1,I2 = sympy.symbols('Vs,R1,R2,Is,Vo,I1,I2')
```

Then we create a list of equations and add one equation for each loop that adds all voltages on the loop


```python
# New list of equations
equations = []

# Loop current equations
equations.append(Vs-R1*I1-R2*(I1-I2)) # Loop current 1
equations.append(-R2*(I2-I1)-Vo)      # Loop current 2
```

Then we create one equation for each current supply that relates it to the loop currents


```python
# Current source equations
equations.append(sympy.Eq(Is,-I2))   # Current source Is
```

Now we can define the unknows for the circuit.

The number of unknowns shall be equal to the number of equations. 

The list includes:

* The loop currents: I1, I2

* The voltage on current sources: Vo

* The measurement values: Vo


```python
# Unknowns list
unknowns = [I1,I2,Vo]
```

We can see the equations and unknows before solving the circuit


```python
showCircuit()
```

    Equations
         -I1*R1 - R2*(I1 - I2) + Vs
         -R2*(-I1 + I2) - Vo
         Eq(Is, -I2)
    
    Unknowns: [I1, I2, Vo]
    


Now we can obtain the solution


```python
# Obtain solution
solution = sympy.solve(equations,unknowns)

print('Vo =',solution[Vo])
```

    Vo = R2*(Is*R1 + Vs)/(R1 + R2)


As in the **nodal** case, you could have used less equations. For instance, you could have used the **Is** current for the second loop.


```python
# Create the circuit symbols
Vs,R1,R2,Is,Vo,I1 = sympy.symbols('Vs,R1,R2,Is,Vo,I1')

# New list of equations
equations = []

# Loop current equations
equations.append(Vs-R1*I1-R2*(I1+Is)) # Loop current 1
equations.append(R2*(Is+I1)-Vo)       # Loop current 2

# Unknowns list
unknowns = [I1,Vo]

# Show equations and unknowns
showCircuit()

# Obtain solution
solution = sympy.solve(equations,unknowns)

print('Vo =',solution[Vo])
```

    Equations
         -I1*R1 - R2*(I1 + Is) + Vs
         R2*(I1 + Is) - Vo
    
    Unknowns: [I1, Vo]
    
    Vo = R2*(Is*R1 + Vs)/(R1 + R2)


## Circuit with current measurements

The following example is special because the circuit has current measurements **I1** and **I2** that we want to obtain



The circuit can be solved using four different methods

### Method #1 : Use nodal method and get currents from resistors

In this method we will just use the normal nodal methods and we will compute the currents using Ohm's law

Note that we don't need the current on **Vs** so there is no point in obtaining the equation on node 1


```python
# Symbols for the circuit
Vs,R1,R2,R3,V2,I1,I2 = sympy.symbols('Vs,R1,R2,R3,V2,I1,I2')

# Nodal equation only on node 2
equations = []
equations.append(-(V2-Vs)/R1-V2/R2-V2/R3)

# Equations for the currents using Ohm's law
equations.append(sympy.Eq(I1,V2/R2))
equations.append(sympy.Eq(I2,V2/R3))

# Unknowns
unknowns = [V2,I1,I2]

# Show equations and unknowns
showCircuit()

# Solve the circuit
solution = sympy.solve(equations,unknowns)

print('I1 =',solution[I1])
print('I2 =',solution[I2])
```

    Equations
         -V2/R3 - V2/R2 + (-V2 + Vs)/R1
         Eq(I1, V2/R2)
         Eq(I2, V2/R3)
    
    Unknowns: [V2, I1, I2]
    
    I1 = R3*Vs/(R1*R2 + R1*R3 + R2*R3)
    I2 = R2*Vs/(R1*R2 + R1*R3 + R2*R3)


### Method #2 : Use four nodes plus ground

We can split the node 2 in three nodes: 2, 3 and 4

That way we can use the current equations to obtain I1 and I2

As all three nodes have the same voltage, we can set them equal using two equations

As in the previous case, we don't need the equation in node 1 


```python
# Symbols for the circuit
Vs,R1,R2,R3,V2,V3,V4,I1,I2 = sympy.symbols('Vs,R1,R2,R3,V2,V3,V4,I1,I2')

# Node equations
equations = []
equations.append(-(V2-Vs)/R1-I1-I2)  # Node 2
equations.append(I1-V3/R2)           # Node 3
equations.append(I2-V4/R3)           # Node 4

# In fact, nodes 2, 3 and 4 are the same
equations.append(sympy.Eq(V2,V3)) 
equations.append(sympy.Eq(V2,V4))

# Unknowns
unknowns = [V2,V3,V4,I1,I2]

# Show equations and unknowns
showCircuit()

# Solve the circuit
solution = sympy.solve(equations,unknowns)

print('I1 =',solution[I1])
print('I2 =',solution[I2])
```

    Equations
         -I1 - I2 + (-V2 + Vs)/R1
         I1 - V3/R2
         I2 - V4/R3
         Eq(V2, V3)
         Eq(V2, V4)
    
    Unknowns: [V2, V3, V4, I1, I2]
    
    I1 = R3*Vs/(R1*R2 + R1*R3 + R2*R3)
    I2 = R2*Vs/(R1*R2 + R1*R3 + R2*R3)


### Method #3 : Use the loop current method

In this case we will define two loop currents 

* Ia goes on the first loop:  Vs -> R1 -> R2
* I2 goes on the second loop: R2 -> R3


```python
# Symbols for the circuit
Vs,R1,R2,R3,Ia,I1,I2 = sympy.symbols('Vs,R1,R2,R3,Ia,I1,I2')

# Loop equations
equations = []
equations.append(Vs-R1*Ia-R2*(Ia-I2))  # Loop Ia
equations.append(-R2*(I2-Ia)-R3*I2)     # Loop I2

# Define I1 from loop currents
equations.append(sympy.Eq(I1,Ia-I2))

# Unknowns
unknowns = [Ia,I1,I2]

# Show equations and unknowns
showCircuit()

# Solve the circuit
solution = sympy.solve(equations,unknowns)

print('I1 =',solution[I1])
print('I2 =',solution[I2])
```

    Equations
         -Ia*R1 - R2*(-I2 + Ia) + Vs
         -I2*R3 - R2*(I2 - Ia)
         Eq(I1, -I2 + Ia)
    
    Unknowns: [Ia, I1, I2]
    
    I1 = R3*Vs/(R1*R2 + R1*R3 + R2*R3)
    I2 = R2*Vs/(R1*R2 + R1*R3 + R2*R3)


### Method #4 : Use a modified loop current method

In this method we will use **I1** and **I2** as loop currents

* I1 goes Vs -> R1 -> R2
* I2 goes Vs -> R1 -> R3


```python
# Symbols for the circuit
Vs,R1,R2,R3,I1,I2 = sympy.symbols('Vs,R1,R2,R3,I1,I2')

# Loop equations
equations = []
equations.append(Vs-R1*(I1+I2)-R2*I1)  # Loop I1
equations.append(Vs-R1*(I1+I2)-R3*I2)  # Loop I2

# Unknowns
unknowns = [I1,I2]

# Show equations and unknowns
showCircuit()

# Solve the circuit
solution = sympy.solve(equations,unknowns)

print('I1 =',solution[I1])
print('I2 =',solution[I2])
```

    Equations
         -I1*R2 - R1*(I1 + I2) + Vs
         -I2*R3 - R1*(I1 + I2) + Vs
    
    Unknowns: [I1, I2]
    
    I1 = R3*Vs/(R1*R2 + R1*R3 + R2*R3)
    I2 = R2*Vs/(R1*R2 + R1*R3 + R2*R3)


## Controlled voltage source

The following circuit includes voltage controlled source



In this case we will treat the controlled voltage source as an independent voltage source, but, we will use $k \cdot V_m$ as its value.

That means that we will need to add an equation to add $V_m$ to the set of symbols

The following code defines and solves the circuit using the loop current method.



```python
# Symbols for the circuit
Vs,R1,R2,k,R3,I1,I2,Vm = sympy.symbols('Vs,R1,R2,k,R3,I1,I2,Vm')

# Loop equations
equations = []
equations.append(Vs-I1*R1-I1*R2)  # Loop I1
equations.append(k*Vm-I2*R3)      # Loop I2

# Equations for Vm and Vo
equations.append(sympy.Eq(Vm,I1*R2))
equations.append(sympy.Eq(Vo,I2*R3))

# Unknowns
unknowns = [I1,I2,Vm,Vo]

# Show equations and unknowns
showCircuit()

# Solve the circuit
solution = sympy.solve(equations,unknowns)

print('Vo =',solution[Vo])
```

    Equations
         -I1*R1 - I1*R2 + Vs
         -I2*R3 + Vm*k
         Eq(Vm, I1*R2)
         Eq(Vo, I2*R3)
    
    Unknowns: [I1, I2, Vm, Vo]
    
    Vo = R2*Vs*k/(R1 + R2)


The circuit could also be solved using the nodal method.

Remember that we don't need the nodal equations for nodes on grounded supplies if we don't need the supply current


```python
# Symbols for the circuit
Vs,Vm,Vo,R1,R2,k = sympy.symbols('Vs,Vm,Vo,R1,R2,k')

# Node equation
equations = []
equations.append(-(Vm-Vs)/R1-Vm/R2)

# Equation for Vo
equations.append(sympy.Eq(Vo,k*Vm))

# Unknowns
unknowns = [Vm,Vo]

# Show equations and unknowns
showCircuit()

# Solve the circuit
solution = sympy.solve(equations,unknowns)

print('Vo =',solution[Vo])
```

    Equations
         -Vm/R2 + (-Vm + Vs)/R1
         Eq(Vo, Vm*k)
    
    Unknowns: [Vm, Vo]
    
    Vo = R2*Vs*k/(R1 + R2)


## RC Circuit

Now we can also solve circuits with capacitor or inductors

You just need to define the **capacitors** currents and voltages as function of the **s** variable:

$\qquad i_C = V_C \cdot C \cdot s \qquad v_C = \frac{i_C}{C \cdot s}$

Also, for **inductors**:

$\qquad i_L = \frac{V_L}{L \cdot s} \qquad v_L = i_L \cdot L \cdot s$

We will use the following circuit as an example



The following code describes and solves the circuit using the current loop method


```python
# Symbols for the circuit
Vs,R1,R3,C1,C2,Vo,I1,I2,s = sympy.symbols('Vs,R1,R3,C1,C2,Vo,I1,I2,s')

# Loop equations
equations = []
equations.append(Vs-R1*I1-(I1-I2)/(C1*s))
equations.append(-(I2-I1)/(C1*s)-R3*I2-I2/(C2*s))

# Equation for Vo
equations.append(sympy.Eq(Vo,I2/(C2*s)))

# Unknowns
unknowns = [I1,I2,Vo]

# Show equations and unknowns
showCircuit()

# Solve the circuit
solution = sympy.solve(equations,unknowns)

Vo_s = solution[Vo]
print('Solution')
print('   Vo =',Vo_s)
```

    Equations
         -I1*R1 + Vs - (I1 - I2)/(C1*s)
         -I2*R3 - I2/(C2*s) + (I1 - I2)/(C1*s)
         Eq(Vo, I2/(C2*s))
    
    Unknowns: [I1, I2, Vo]
    
    Solution
       Vo = -C1*Vs/(C2 - (C1*R1*s + 1)*(C1*C2*R3*s + C1 + C2))


We can obtain a better solution equation using **simplify**, **expand** and **collect** from the sympy module


```python
# Use simplify, expand and collect to give a prettier equation
Vo_s = Vo_s.simplify().expand()          # Eliminate quotients of quotients
Vo_s = sympy.collect(Vo_s,s).simplify()  # Group s symbols and simplify

print('Prettier solution')
print('   Vo =',Vo_s)
```

    Prettier solution
       Vo = Vs/(C1*C2*R1*R3*s**2 + s*(C1*R1 + C2*R1 + C2*R3) + 1)


We can obtian a particular solution using numbers to substitute the literals

To substitute one symbol you can use:

>`expr.subs(oldSymbol,newSymbol)`

But, in order to substitute several symbols at once you can use a substitution dictionary:

>`expr.subs({old1:new1,old2:new2,....})`

Substituting in our example we get a **particular** solution


```python
H_s = Vo_s.subs({Vs:1,R1:1000,R3:100,C1:1e-6,C2:100e-9})
H_s = H_s.simplify()

print('Particular solution')
print('    H(s) = Vo(s)/Vs(s) =',H_s)
```

    Particular solution
        H(s) = Vo(s)/Vs(s) = 1/(1.0e-8*s**2 + 0.00111*s + 1)


We can also get the **poles**, the **zeros** and the **DC gain**


```python
numer,denom =H_s.as_numer_denom()
print('Num =',numer)
print('Den =',denom)
print()

zeros = sympy.roots(numer,s)
poles = sympy.roots(denom,s)
print('Zeros =',zeros)
print('Poles =',poles)
print()

print('DC gain =',H_s.subs(s,0).evalf())
```

    Num = 1
    Den = 1.0e-8*s**2 + 0.00111*s + 1
    
    Zeros = {}
    Poles = {-110091.666030631: 1, -908.333969368548: 1}
    
    DC gain = 1.00000000000000


## Opamp circuit

We can also solve operational amplifier circuits



The easiest solution is obtained if we can guarantee that the **virtual shortcircuit** holds

In this case, the opamp output is an unknown and the two input voltages are equal:

$\qquad V_{(+)}=V_{(-)}$ 

We will use the nodal method.

Remember that we don't need the node equations in nodes 1 and 3 because we don't need the source currents


```python
# Symbols for the circuit
Vs,Ri,Rf,Vo,V2 = sympy.symbols('Vs,Ri,Rf,Vo,V2')

# Node equation
equations = []
equations.append(-(V2-Vs)/Ri-(V2-Vo)/Rf)

# Virtual shortcircuit
equations.append(V2)  # V2 = V(+) = 0

# Unknowns
unknowns = [V2,Vo]

# Show equations and unknowns
showCircuit()

# Solve the circuit
solution = sympy.solve(equations,unknowns)

print('Vo =',solution[Vo])
```

    Equations
         (-V2 + Vs)/Ri - (V2 - Vo)/Rf
         V2
    
    Unknowns: [V2, Vo]
    
    Vo = -Rf*Vs/Ri


### Finite gain solution

If we don't want to use the **virtual short circuit** we can solve the opamp as a voltage controlled voltage source

In this case the opamp can be defined with two equations:

$\qquad V_d = V_{(+)}-V_{(-)}$

$\qquad  V_O = A \cdot V_d $

Then, the ideal case will be given when $A \rightarrow \infty$


```python
# Symbols for the circuit
Vs,Ri,Rf,Vo,V2,Vd,A = sympy.symbols('Vs,Ri,Rf,Vo,V2,Vd,A')

# Node equation
equations = []
equations.append(-(V2-Vs)/Ri-(V2-Vo)/Rf)

# Opamp equations
equations.append(sympy.Eq(Vd,0-V2))
equations.append(sympy.Eq(Vo,A*Vd))

# Unknowns
unknowns = [V2,Vo,Vd]

# Show equations and unknowns
showCircuit()

# Solve the circuit
solution = sympy.solve(equations,unknowns)

print('Solution as function of A')
print()
print('   Vo =',solution[Vo])
print()
print('Solution for A -> oo')
print()
print('   Vo =',sympy.limit(solution[Vo],A,sympy.oo))
```

    Equations
         (-V2 + Vs)/Ri - (V2 - Vo)/Rf
         Eq(Vd, -V2)
         Eq(Vo, A*Vd)
    
    Unknowns: [V2, Vo, Vd]
    
    Solution as function of A
    
       Vo = -A*Rf*Vs/(A*Ri + Rf + Ri)
    
    Solution for A -> oo
    
       Vo = -Rf*Vs/Ri


### Dominant pole solution

Having a solution as function of **A** enables us to obtain the response of the circuit using a dominant pole model for the operational amplifier.

Just susbtitute, in the solution **A** for the one pole model of the opamp

$\qquad A = \frac{Ao \cdot p1}{s+p1}$


```python
# New A, p1 and s symbols
Ao,p1,s = sympy.symbols('Ao,p1,s')

Vo_s = solution[Vo].subs(A,Ao*p1/(s+p1))

# Use simplify, expand and collect to give a prettier equation
Vo_s = Vo_s.simplify().expand()  # Eliminate quotients of quotients
Vo_s = sympy.collect(Vo_s,s)     # Group s symbols
Vo_s = sympy.collect(Vo_s,Ri)    # Group Ri symbols

print('Vo(s) =',Vo_s)
```

    Vo(s) = -Ao*Rf*Vs*p1/(Rf*p1 + Ri*(Ao*p1 + p1) + s*(Rf + Ri))


We can obtian, as in a previous example, a particular solution using numbers to substitute the literals

In our opamp circuit solution we substitute:

* Circuit resistors: R1, R2

* Opamp model: Ao, p1

We also substitute Vs for 1 to obtain the transfer function $H(s)$

$\qquad H(s)=\frac{V_O(s)}{V_s(s)}$


```python
H_s = Vo_s.subs({Vs:1,Ao:100000,Rf:100000,Ri:10000,p1:16})
print('H(s) =',H_s)
```

    H(s) = -160000000000/(110000*s + 16001760000)


Now you can also obtain the **poles** and **zeros** of $H(s)$

Also we can get the DC gain

We will use the **evalf()** method that evaluates a **sympy** expression to a floating point number


```python
numer,denom =H_s.as_numer_denom()
print('Num =',numer)
print('Den =',denom)
print()

zeros = sympy.roots(numer,s)
poles = sympy.roots(denom,s)
print('Zeros =',zeros)
print('Poles =',poles)
print()

print('DC gain =',H_s.subs(s,0).evalf())
```

    Num = -160000000000
    Den = 110000*s + 16001760000
    
    Zeros = {}
    Poles = {-1600176/11: 1}
    
    DC gain = -9.99890012098669


<BR><BR><BR><BR><BR><BR>

## Document information

Copyright © Vicente Jiménez (2018)

Last update: 26/3/2018

This work is licensed under a [Creative Common Attribution-ShareAlike 4.0 International license](http://creativecommons.org/licenses/by-sa/4.0/). 

You can find the module [here](https://github.com/R6500/Python-bits/tree/master/Modules)

See my blogs [AIM65](http://aim65.blogspot.com.es/) (in spanish) and [R6500](http://r6500.blogspot.com.es/) (in english)


```python

```

9abdfdb1a23a795342af855c1210a656867e1b1a
ipynb
Jupyter Notebook
notebooks/sympy-circuit-sandbox.ipynb
IronwoodLabs/ambrose
b5f3570b8a61ad22cbd9bf95b7217bcec118e8f4
MIT
notebooks/sympy-circuit-sandbox.ipynb
IronwoodLabs/ambrose
b5f3570b8a61ad22cbd9bf95b7217bcec118e8f4
MIT
notebooks/sympy-circuit-sandbox.ipynb
IronwoodLabs/ambrose
b5f3570b8a61ad22cbd9bf95b7217bcec118e8f4
MIT
Qwen/Qwen-72B
1. YES
2. YES