Challenging RLC Circuit

The state variables are the capacitor voltages and the inductor current $(V_{C1}, V_{C2}, I_L)$ . Put together a system of equations which relates the time-derivates of the state variables to the state variables and the forcing function. We need three equations to solve for the time-derivatives of the state variables. You can write the equations however you like, as long as the system matrix ends up being invertible. The first equation is a KVL loop around the outer edge (the second line is an expanded version of the first line).

$V_S - R_1 I_{C1} - V_{C1} - R_4 I_{C2} - V_{C2} = 0 \\ V_S - R_1 C_1 \dot{V}_{C1} - V_{C1} - R_4 C_2 \dot{V}_{C2} - V_{C2} = 0$

The second equation relates the currents flowing out of the source to the currents flowing back into the source (the second line is an expanded version of the first line).

$I_L + I_{C1} = I_{R3} + I_{C2} \\ I_L + C_1 \dot{V}_{C1} = \frac{V_S - L \dot{I}_L}{R3} + C_2 \dot{V}_{C2}$

The third equation relates two different expressions for the voltage across $R_2$ .

$(V_S - R_1 C_1 \dot{V}_{C1} - V_{C1}) - (V_S - L \dot{I}_L) = R_2 I_{R2} = R_2 (C_1 \dot{V}_{C1} - C_2 \dot{V}_{C2})$

Re-arranging these equations somewhat results in the following:

$\begin{pmatrix} R_1 C_1 & R_4 C_2 & 0 \\ R_3 C_1 & -R_3 C_2 & L \\ R_2 C_1 + R_1 C_1 & -R_2 C_2 & -L \\ \end{pmatrix} \begin{pmatrix} \dot{V}_{C1} \\ \dot{V}_{C2} \\ \dot{I}_L \\ \end{pmatrix} = \begin{pmatrix} V_S - V_{C1} - V_{C2} \\ V_S - R_3 I_L \\ -V_{C1} \\ \end{pmatrix}$

Solve the linear system on each time step for the time-derivatives of the state variables, and numerically integrate. The source current is:

$I_S = I_L + C_1 \dot{V}_{C1}$

A plot of the source current is shown below. The current starts at $60$ initially and asymptotically approaches $100$ . It has a minimum value of about $6.85$ . I also double-checked this result using the Electromagnetic Transients Program (EMTP).

Python simulation code:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94

import math
import numpy as np

VS = 10.0

R1 = 0.1
R2 = 0.1
R3 = 0.1
R4 = 0.1

L = 1.0

C1 = 1.0
C2 = 2.0

##################################

J11 = R1*C1
J12 = R4*C2
J13 = 0.0

J21 = R3*C1
J22 = -R3*C2
J23 = L

J31 = R2*C1 + R1*C1
J32 = -R2*C2 
J33 = -L 

J = np.array([[J11,J12,J13],[J21,J22,J23],[J31,J32,J33]])

##################################

dt = 10.0**(-4.0)
t = 0.0
count = 0

IL = 0.0
VC1 = 0.0
VC2 = 0.0

###########################################

vec = np.array([VS-VC1-VC2,VS-R3*IL,-VC1])

Sol = np.linalg.solve(J,vec)

VC1d = Sol[0]
VC2d = Sol[1]
ILd = Sol[2]

IS = IL + C1*VC1d

#print t,IS

###########################################

#print "t VS Vout"

ISmin = 99999999.0

while t <= 100.0:

    VC1 = VC1 + VC1d * dt
    VC2 = VC2 + VC2d * dt
    IL = IL + ILd * dt

    vec = np.array([VS-VC1-VC2,VS-R3*IL,-VC1])

    Sol = np.linalg.solve(J,vec)

    VC1d = Sol[0]
    VC2d = Sol[1]
    ILd = Sol[2]

    IS = IL + C1*VC1d

    t = t + dt
    count = count + 1

    #if count % 1000 == 0:
        #print t,IS

    if IS < ISmin:
        ISmin = IS

##################################

print ""
print ""

print dt
print ISmin
print (160.0/ISmin)