RLC (4-24-2021)

I will post my solution as well, which is the same as that of @Karan Chatrath, apart from minor stylistic differences. The key here is to write the time derivatives of the state variables in terms of the state variables and the forcing function (that is known as a state-space formulation). The state variables are the inductor currents and the capacitor voltages, and the forcing function is the DC source. Calculating the time derivatives allows for numerical integration.

  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
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143

# Constants

dt = 10.0**(-5.0)

VS = 10.0

R1 = 1.0
R2 = 1.0
R3 = 1.0

C1 = 5.0
C2 = 5.0

L1 = 1.0
L2 = 1.0
L3 = 1.0

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

# Initialize state variables

t = 0.0
count = 0

IL1 = 0.0
IL2 = 0.0
IL3 = 0.0

VC1 = 0.0
VC2 = 0.0

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

# Initialize time derivatives and source current

VL1 = VS - VC2 - VC1
IL1d = VL1/L1  #

IR1 = (VS - VL1)/R1

IC1 = IL1 - IR1

VC1d = IC1/C1  #

VR2 = VS - VC2
IR2 = VR2/R2

IC2 = IC1 + IR2 + IL2
VC2d = IC2/C2  #

IR3 = IL2 + IL3
VR3 = R3*IR3

VL2 = VR2 - VR3
IL2d = VL2/L2  #

VL3 = VS - VR3
IL3d = VL3/L3  #

IS = IR1 + IC2 + IL3

#print t,IS

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

ISmax = 0.0

# Time simulation

while t <= 100.0:

    IL1 = IL1 + IL1d*dt     # Explicit Euler integration
    IL2 = IL2 + IL2d*dt
    IL3 = IL3 + IL3d*dt

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

    ########

    VL1 = VS - VC2 - VC1   # Calculate time derivatives and source current
    IL1d = VL1/L1  #        # State space

    IR1 = (VS - VL1)/R1

    IC1 = IL1 - IR1

    VC1d = IC1/C1  #

    VR2 = VS - VC2
    IR2 = VR2/R2

    IC2 = IC1 + IR2 + IL2
    VC2d = IC2/C2  #

    IR3 = IL2 + IL3
    VR3 = R3*IR3

    VL2 = VR2 - VR3
    IL2d = VL2/L2  #

    VL3 = VS - VR3
    IL3d = VL3/L3  #

    IS = IR1 + IC2 + IL3

    ########

    if IS > ISmax:
        ISmax = IS

    t = t + dt
    count = count + 1

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

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

print ""
print ""

print dt
print t
print ""
print ISmax
print (ISmax/(10.0 + 30.0))

#>>> 
#0.0001
#100.000000002

#36.3235594067
#0.908088985168

#>>> ================================ RESTART ================================
#>>> 
#1e-05
#100.000000021

#36.3234065033
#0.908085162582
#>>> 

The circuit equations as per Kirchoff's current and voltage laws are:

$I_{C1} +I_{R1} = I_{L1} \ \dots(1)$ $I_{C2}= I_{C1} + I_{L2} + I_{R2} \ \dots(2)$ $I_{R3}=I_{L3} + I_{L2} \ \dots(3)$ $I_S = I_{R1}+I_{C2}+I_{L3} \ \dots(4)$ $\dot{I}_{L1} = V_S - I_{R1} \ \dots(5)$ $\dot{Q}_{C1} = I_{C1} \ \dots(6)$ $\dot{Q}_{C2} = I_{C2} \ \dots(7)$ $\dot{I}_{L2} = I_{R2} - I_{R3} \ \dots(8)$ $I_{R1} = \frac{Q_{C1} + Q_{C2}}{5} \ \dots(9)$ $Q_{C1} + \dot{I}_{L1} = I_{R2} \ \dots(10)$ $\dot{I}_{L3} - \dot{I}_{L2} - \frac{Q_{C2}}{5}=0 \ \dots(11)$

Initial conditions:

$I_{L1}(0) =I_{L2}(0)=I_{L3}(0) = Q_{C1}(0)= Q_{C2}(0) = 0$

By inspecting the circuit, one can easily conclude that at $t=0$ the source current $I_S = 10 \ A$ and at steady state $I_S = 30 \ A$ . These results are found using the initial conditions.

A plot of source current as a function of time is as follows:

This is a 5th order circuit. Instead of simplifying the expressions, I have chosen a purely numerical approach:

 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

clear all
clc

% Initial conditions:
Q_C1(1)     = 0;
Q_C2(1)     = 0;
I_L1(1)     = 0;
I_L2(1)     = 0;
I_L3(1)     = 0;

% Source Voltage:
VS          = 10;

% Time step, final time and time array:
dt          = 1e-4;
tf          = 100;
t           = 0:dt:tf;

% Loop to simulate the circuit:
for k = 1:length(t)-1

  I_R1          = (Q_C1(k) + Q_C2(k))/5;      % Equation 9

  dI_L1         = VS - I_R1;                  % Equation 5
  I_L1(k+1)     = I_L1(k) + dI_L1*dt;         % Explicit Euler

  I_R2          = dI_L1 + Q_C1(k)/5;          % Equation 10

  I_C1          = I_L1(k) - I_R1;             % Equation 1
  dQ_C1         = I_C1;                       % Equation 6
  Q_C1(k+1)     = Q_C1(k) + dt*dQ_C1;         % Explicit Euler

  I_C2          = I_C1 + I_R2 +I_L2(k);       % Equation 2
  dQ_C2         = I_C2;                       % Equation 7
  Q_C2(k+1)     = Q_C2(k) + dt*dQ_C2;         % Explicit Euler

  I_R3          = I_L2(k) + I_L3(k);          % Equation 3
  dI_L2         = I_R2 - I_R3;                % Equation 8
  dI_L3         = dI_L2 + Q_C2(k)/5;          % Equation 11

  I_L2(k+1)     = I_L2(k) + dt*dI_L2;         % Explicit Euler
  I_L3(k+1)     = I_L3(k) + dt*dI_L3;         % Explicit Euler

  I_S(k)        = I_R1 + I_C2 + I_L3(k);      % Equation 4

  % Computing source current at final time step to prevent array length mismatches:
  if k == length(t)-1
    I_S(k+1)    = I_R1 + I_C2 + I_L3(k+1);    % Equation 4
  end

end

% Required Answer:
ANSWER   = max(I_S)/(10 + 30)
%ANSWER  ~ 0.9081