Question 1

Determine the S a2 duty ratio in order to the dc bus voltage would be 400 volts, while the fuel
cell voltage is equal to its nominal value and the power necessary is its rated value (Consider
that the inductance is high enough in such a way that the input current has no ripple.).

In [1]: from sympy import *

import numpy as np
import [Link] as plt

First way
The easiest way to solve this question is to consider that the average voltage value over the
indutor during the switching period is zero in stead state and the output capacitor voltage is
ripple free or the output capacitance is very enough to desconsider any ripple. From these
assumptions, it is possible to state the following equation:

Ts
1
⟨vL ⟩ = ∫ vL dt = 0
Ts 0

The voltage during the first operation state, S closed, is equal to the V
2 in
and during the
second operation state, S closed , is equal to the (V
1 in − Vout ) .

DTs Ts
1
(∫ Vin dt + ∫ (Vin − Vout )dt) = 0
Ts
0 DTs

DTs Vin + (Ts − DTs )(Vin − Vout ) = 0

Vout 1
=
Vin 1 − D

Considering the input voltage equal to 200 volts and the desired output voltage to be 400 volts,
the above equation becomes :

400V 1 1
= 2 = → 1 − D =
200V 1 − D 2

And the duty cycle value is:

D = 0.5
Another way to find the duty cycle is by solving a linear equation system considering the
specification of all components.

Second Way
This method consists of calculating the average value of the output capacitor directly. For this, it
is necessary to know what are the equation that describes the voltage waveform first and then
calculate its average value over a switching period.

Waveform equation
The assumptions considered for this method are:

The components are linear;

The capacitor voltage and the inductor current in the end of one topological state have the
same value at the beginning of the following topological state.

There are two topological states in this case:

One for the S transistor closed while S is open;

2 1

and one for the S transistor closed while S is open;

1 2

For each state, there is one equivalent circuit and two state equations.

t
A1 t A1 (t−τ )
x1 (t) = e x1 (0) + ∫ e B1 u(t)dτ
0

t
A2 t A2 (t−τ )
x2 (t) = e x2 (0) + ∫ e B2 u(t)dτ
0

The state values at the end of each operation stage is given by:

D 1 Ts
A 1 D 1 Ts A1 (t−τ )
x1 (D1 Ts ) = e x1 (0) + ∫ e B1 u(t)dτ
0

D 2 Ts
A 2 D 2 Ts A2 (t−τ )
x2 (D2 Ts ) = e x2 (0) + ∫ e B2 u(t)dτ
0

From the second assumption, these relations can be made:

x2 (0) = x1 (D1 Ts )

x1 (0) = x2 (D2 Ts )

And the system can be rewritten as:

D 1 Ts
A 1 D 1 Ts A1 (t−τ )
x2 (0) − e x1 (0) = ∫ e B1 u(t)dt
0

A 1 D 1 Ts D 2 Ts A2 (t−τ )
x1 (0) − e x2 (0) = ∫ e B2 u(t)dt
0
 or as in the matricial form:

D 1 Ts
A1 (t−τ )
−e
A 1 D 1 Ts
I x1 (0) ⎡∫ e B1 u(τ )dτ ⎤
0
[ ][ ] =
A 2 D 2 Ts D 2 Ts
I −e x2 (0) ⎣ A2 (t−τ ) ⎦
∫ e B2 u(τ )dτ
0

State systems
In this case, the components specifications influence the result. The following values will be
used:

In [2]: L = 100e-6
C = 4e-6
R = 400**2/75e3
E = 200
Ts = 10e-6

In [3]: t = symbols('t', positive = true)

tau = symbols('tau', positive = true)

D1 = symbols('D1', positive = true)

D2 = symbols('D2', positive = true)

The next step is the determination of matrix A and vector B that describe the equivalent circuit
for each operational state.

dx
= Ax + Bu(t)
dt

First State System

From the circuit, the following equation can be written:

Inductor current equation

dil dil E
vl = E → = E → =
dt dt L

Capacitor voltage equation

vc dvc vc dvc vc
ic = − → C = − → = −
R dt R dt RC

The subscription 1 refers to the first operation stage.

dx1
= A1 x + B1 u1 (t)
dt

il
x1 = [ ]
vc
 1
dx1 0 0
L
= [ 1
] x1 + [ ][E]
dt 0 −
RC
0

Matrix A

In [4]: A1 = Matrix([[0,0],
[0,-1/(R*C)]])
A1

Out[4]: 0 0
[ ]
0 −117187.5

Vector B

In [5]: B1 = Matrix([[1/L],
[0]])
B1

Out[5]: 10000.0
[ ]
0

Input

In [6]: u1 = E

Second State System

From the circuit, the following equation can be written:

Inductor current equation

dil dil 1 E
vl = E − vc → = E − vc → = − vc +
dt dt L L

Capacitor voltage equation

vc dvc vc dvc il vc
ic = il − → C = il − → = −
R dt R dt C RC

dx2
= A2 x2 + B2 u(t)
dt

The subscription 2 refers to the second operation stage.

il
x2 = [ ]
vc

1
1
0 −
dx2 L
L
= [ ] x1 + [ ][E]
1 1
dt − 0
C RC
 Matrix A

In [7]: A2 = Matrix([[0,-1/L],
[1/C,-1/(R*C)]])
A2

Out[7]: 0 −10000.0
[ ]
250000.0 −117187.5

Vector B

In [8]: B2 = Matrix([[1/L],
[0]])
B2

Out[8]: 10000.0
[ ]
0

Input

In [9]: u2 = E

Equation System Ax = b

D 1 Ts
A1 (t−τ )
−e
A 1 D 1 Ts
I x1 (0) ⎡∫ e B1 u(τ )dτ ⎤
0
[ ][ ] =
A 2 D 2 Ts D 2 Ts
I −e x2 (0) ⎣ A2 (t−τ ) ⎦
∫ e B2 u(τ )dτ
0

The next step is to calculate the exponential matrices in the system below:

Matrix e A1 t

In [10]: expA1t = exp(A1*(t)).expand(complex=True)

expA1t

Out[10]: 1.0 0
[ ]
−117187.5t
0 1.0e

Matrix e A2 t

In [11]: expA2t = exp(A2*(t)).expand(complex=True)

N(expA2t,3)

Out[11]: −0.459e
−89142.5231187768t
+ 1.46e
−28044.9768812232t
0.164e
−89142.5231187768t
− 0.164e
−28044.97688

[
−89142.5231187768t −28044.9768812232t −89142.5231187768t −28044.976881
−4.09e + 4.09e 1.46e − 0.459e

Vector ∫
D 1 Ts
A1 (t−τ )
e B1 u(τ )dτ
0
In [12]: expA1ttau = exp(A1*(t-tau)).expand(complex=True)
expA1ttau

Out[12]: 1.0 0
[ ]
−117187.5t 117187.5τ
0 1.0e e

In [13]: int_expA1B1u1_0_D1Ts = integrate((expA1ttau*B1*u1).expand(complex=True),(tau,0,t)).sub

int_expA1B1u1_0_D1Ts

Out[13]: 20.0D1
[ ]
0

Vector ∫
D 2 Ts
A2 (t−τ )
e B2 u(τ )dτ
0

In [14]: expA2ttau = exp(A2*(t-tau)).expand(complex=True)

N(expA2ttau,4)

Out[14]: −0.459e
−89142.5231187768t
e
89142.5231187768τ
+ 1.459e
−28044.9768812232t
e
28044.9768812232τ
0.1637e
−8

[
−89142.5231187768t 89142.5231187768τ −28044.9768812232t 28044.9768812232τ −8
−4.092e e + 4.092e e 1.459e

In [15]: int_expA2B2u2_0_D1Ts = integrate((expA2ttau*B2*u2).expand(complex=True),(tau,0,t)).sub

N(int_expA2B2u2_0_D1Ts,4)

Out[15]: −78.6e
0.280449768812232D1
+ 4.223e
0.891425231187768D1
+ 93.75
[ ]
0.280449768812232D1 0.891425231187768D1
−220.4e + 37.65e + 200.0

A matrix

In [16]: A = Matrix([[-[Link](t,D1*Ts),eye(2)],[eye(2),-[Link](t,(1-D1)*Ts)]])
N(A,4)

Out[16]: ⎡
−1.0 0 1.0

−1.171875D1
⎢ 0 −1.0e 0
⎢
⎢ 0.280449768812232D1 0.891425231187768D1 0.280
⎢ 1.0 0 −1.102e + 0.1882e 0.1236e

⎣ 0.280449768812232D1 0.891425231187768D1 0.28

0 1.0 −3.091e + 1.678e 0.3468e

B Matrix

In [17]: B = Matrix([int_expA1B1u1_0_D1Ts ,int_expA2B2u2_0_D1Ts ] )

N(B,4)

Out[17]: ⎡
20.0D1
⎤

⎢ 0 ⎥
⎢ ⎥
⎢ 0.280449768812232D1 0.891425231187768D1 ⎥
⎢ −78.6e + 4.223e + 93.75 ⎥

⎣ 0.280449768812232D1 0.891425231187768D1 ⎦
−220.4e + 37.65e + 200.0

Initial Values
In [18]: Xinitial = [Link](D1,0.5).inv()*[Link](D1,0.5)
N(Xinitial,4)

Out[18]: ⎡
358.5
⎤

⎢ 498.8 ⎥
⎢ ⎥
⎢ 368.5 ⎥

⎣ ⎦
277.6

In [19]: X10 = Matrix([[1,0,0,0],[0,1,0,0]])*Xinitial

N(X10,4)

Out[19]: 358.5
[ ]
498.8

In [20]: X20 = Matrix([[0,0,1,0],[0,0,0,1]])*Xinitial

N(X20,4)

Out[20]: 368.5
[ ]
277.6

Waveforms equations for the first operation stage

In [21]: Y1 = expA1t*X10 + integrate((expA1ttau*B1*u1).expand(complex=True),(tau,0,t))

N(Y1,3)

Out[21]: 6
2.0 ⋅ 10 t + 359.0
[ ]
−117187.5t
499.0e

Waveforms equations for the second operation stage

In [22]: Y2 = expA2t*X20 + integrate((expA2ttau*B2*u2).expand(complex=True),(tau,0,t))

N(Y2,3)

Out[22]: 93.8 − 113.0e

−89142.5231187768t
+ 388.0e
−28044.9768812232t

[ ]
3 −89142.5231187768t 3 −28044.9768812232t
200.0 − 1.01 ⋅ 10 e + 1.09 ⋅ 10 e

Average Value
In [23]: average_values = integrate(Y1,(t,0,0.5*Ts))/Ts + integrate(Y2,(t,0,0.5*Ts))/Ts
N(average_values,4)

Out[23]: 364.0
[ ]
388.7

The actual average voltage value for a duty cycle of 0.5 is 388 volts and not 400 volts as
required, at least for the chosen components specifications. This difference is due to the
 assumption that the capacitor voltage waveform was considerated without ripple in the first
method while its real waveform exponencial behavior was considered here.

Below, the same method is implemented as a function and results are obtained for the same
components specifications values, except for the output capacitance.

In [24]: def fun(D,*args):

R,L,C,E,Ts = args

A1 = Matrix([[0,0],[0,-1/(R*C)]])
B1 = Matrix([[1/L],[0]])
u1 = E

A2 = Matrix([[0,-1/L],[1/C,-1/(R*C)]])
B2 = Matrix([[1/L],[0]])
u2 = E

expA1t = exp(A1*(t)).expand(complex=True)
expA2t = exp(A2*(t)).expand(complex=True)

expA1ttau = exp(A1*(t-tau)).expand(complex=True)
expA2ttau = exp(A2*(t-tau)).expand(complex=True)

int_expA1B1u1_0_D1Ts = integrate((expA1ttau*B1*u1).expand(complex=True),(tau,0,t))
int_expA2B2u2_0_D2Ts = integrate((expA2ttau*B2*u2).expand(complex=True),(tau,0,t))

A = Matrix([[-[Link](t,D*Ts),eye(2)],[eye(2),-[Link](t,(1-D)*Ts)]])
B = Matrix([int_expA1B1u1_0_D1Ts ,int_expA2B2u2_0_D2Ts ] )

Xinitial = [Link](D,0.5).inv()*[Link](D,0.5)

X10 = Matrix([[1,0,0,0],[0,1,0,0]])*Xinitial
X20 = Matrix([[0,0,1,0],[0,0,0,1]])*Xinitial

Y1 = expA1t*X10 + integrate((expA1ttau*B1*u1).expand(complex=True),(tau,0,t))
Y2 = expA2t*X20 + integrate((expA2ttau*B2*u2).expand(complex=True),(tau,0,t))

average_values = integrate(Y1,(t,0,D*Ts))/Ts + integrate(Y2,(t,0,(1-D)*Ts))/Ts

return average_values

The results are obtained for 4, 8 and 40 μF.

In [25]: N(fun(0.5,*(400**2/75e3,100e-6,4e-6,200,10e-6 )),4)

Out[25]: 364.0
[ ]
388.7

In [26]: N(fun(0.5,*(400**2/75e3,100e-6,8e-6,200,10e-6 )),4)

Out[26]: 371.9
[ ]
396.9
In [27]: N(fun(0.5,*(400**2/75e3,100e-6,16e-6,200,10e-6 )),4)

Out[27]: 374.1
[ ]
399.2

As can be noted, as the capacitance value increases, the average voltage value gets close to 400
volts.

In [28]: D_array = [Link](0.49,0.51,10)

Vc_4uF_array = [fun(d,*(400**2/75e3,100e-6,4e-6,200,10e-6 ))[1] for d in D_array]

Vc_8uF_array = [fun(d,*(400**2/75e3,100e-6,8e-6,200,10e-6 ))[1] for d in D_array]
Vc_40uF_array = [fun(d,*(400**2/75e3,100e-6,40e-6,200,10e-6 ))[1] for d in D_array]

Vc_ideal_array = [200*1/(1-d) for d in D_array]

The graph below presents the average capacitor voltage value around 50% duty cycle for 4
cases:

capacitance of 4μF
capacitance of 8μF
capacitance of 40μF
Ideal or infinite capacitance.

In [29]: [Link](D_array,Vc_4uF_array,label='C=$4μF/ΔV%={{:3.4}}$'.format( (1-exp(-5e-6/(R*4e-

[Link](D_array,Vc_8uF_array,label='C=$8μF/ΔV%={:3.4}$'.format( (1-exp(-5e-6/(R*8e-6)
[Link](D_array,Vc_40uF_array,label='C=$40μF/ΔV%={:3.4}$'.format( (1-exp(-5e-6/(R*40e
[Link](D_array,Vc_ideal_array,'k-',label='Ideal')

[Link]()
[Link]()

[Link]('$D$')
[Link]('$\\left < v_c \\right >$ [V]')

Out[29]: Text(0, 0.5, '$\\left < v_c \\right >$ [V]')

 If the capacitance is high enough, the actual value and the one calculated by the first method
are very similar.

In [30]: def fun_graph(D,*args):

R,L,C,E,Ts = args

A1 = Matrix([[0,0],[0,-1/(R*C)]])
B1 = Matrix([[1/L],[0]])
u1 = E

A2 = Matrix([[0,-1/L],[1/C,-1/(R*C)]])
B2 = Matrix([[1/L],[0]])
u2 = E

expA1t = exp(A1*(t)).expand(complex=True)
expA2t = exp(A2*(t)).expand(complex=True)

expA1ttau = exp(A1*(t-tau)).expand(complex=True)
expA2ttau = exp(A2*(t-tau)).expand(complex=True)

int_expA1B1u1_0_D1Ts = integrate((expA1ttau*B1*u1).expand(complex=True),(tau,0,t))
int_expA2B2u2_0_D2Ts = integrate((expA2ttau*B2*u2).expand(complex=True),(tau,0,t))

A = Matrix([[-[Link](t,D*Ts),eye(2)],[eye(2),-[Link](t,(1-D)*Ts)]])
B = Matrix([int_expA1B1u1_0_D1Ts ,int_expA2B2u2_0_D2Ts ] )

Xinitial = [Link](D,0.5).inv()*[Link](D,0.5)

X10 = Matrix([[1,0,0,0],[0,1,0,0]])*Xinitial
 X20 = Matrix([[0,0,1,0],[0,0,0,1]])*Xinitial

Y1 = expA1t*X10 + integrate((expA1ttau*B1*u1).expand(complex=True),(tau,0,t))
Y2 = expA2t*X20 + integrate((expA2ttau*B2*u2).expand(complex=True),(tau,0,t))

t_array = [Link](0,Ts,100)

vc_array = [Y1[1].subs(t,tx) if tx < D*Ts else Y2[1].subs(t,tx-D*Ts) for tx in t_a

il_array = [Y1[0].subs(t,tx) if tx < D*Ts else Y2[0].subs(t,tx-D*Ts) for tx in t_a

return t_array,il_array,vc_array

Below, the waveforms for the three finite capacitance values are plotted. It can be observed that
as the capacitance increases, the ripple decreases and the voltage waveform slope becomes
more constant over the operation stage interval. This behavior let the average voltage value to
tend to the value calculate in the first method.

In [31]: t_array,il_array_4uF,vc_array_4uF = fun_graph(0.5,*(400**2/75e3,100e-6,4e-6,200,10e-6

t_array,il_array_8uF,vc_array_8uF = fun_graph(0.5,*(400**2/75e3,100e-6,8e-6,200,10e-6
t_array,il_array_40uF,vc_array_40uF = fun_graph(0.5,*(400**2/75e3,100e-6,40e-6,200,10e

In [32]: [Link](t_array*1e6,vc_array_4uF,label='C=$4μF $')

[Link](t_array*1e6,vc_array_8uF,label='C=$8μF $')
[Link](t_array*1e6,vc_array_40uF,label='C=$40μF $')

[Link]('$t \, [μs]$')
[Link]('$v_c$ [V]')

[Link]()

[Link](0,10)
[Link](0,600)
[Link]()
 The output filter is designed to have a very low ripple, usually 1% or less of the nominal voltage,
hence, the first method voltage gain equation is enough to have a very good idea of the actual
duty cycle value. Also, components non-idealities have a higher impact over the duty cycle.

For example, the function fun was altered to consider a inductor resistance of 0.5% of the
output load, (10.7mΩ). This is sufficient to decrease the average voltage from 399.2 to 391.4
volts.

In [33]: def fun_r(D,*args):

R,L,C,E,Ts = args
rl=0.005*R

A1 = Matrix([[-rl/L,0],[0,-1/(R*C)]])
B1 = Matrix([[1/L],[0]])
u1 = E

A2 = Matrix([[-rl/L,-1/L],[1/C,-1/(R*C)]])
B2 = Matrix([[1/L],[0]])
u2 = E

expA1t = exp(A1*(t)).expand(complex=True)
expA2t = exp(A2*(t)).expand(complex=True)

expA1ttau = exp(A1*(t-tau)).expand(complex=True)
expA2ttau = exp(A2*(t-tau)).expand(complex=True)

int_expA1B1u1_0_D1Ts = integrate((expA1ttau*B1*u1).expand(complex=True),(tau,0,t))
 int_expA2B2u2_0_D2Ts = integrate((expA2ttau*B2*u2).expand(complex=True),(tau,0,t))

A = Matrix([[-[Link](t,D*Ts),eye(2)],
[eye(2),-[Link](t,(1-D)*Ts)]])

B = Matrix([int_expA1B1u1_0_D1Ts ,
int_expA2B2u2_0_D2Ts ] )

Xinitial = [Link](D,0.5).inv()*[Link](D,0.5)

X10 = Matrix([[1,0,0,0],
[0,1,0,0]])*Xinitial

X20 = Matrix([[0,0,1,0],
[0,0,0,1]])*Xinitial

Y1 = expA1t*X10 + integrate((expA1ttau*B1*u1).expand(complex=True),(tau,0,t))
Y2 = expA2t*X20 + integrate((expA2ttau*B2*u2).expand(complex=True),(tau,0,t))

average_values = integrate(Y1,(t,0,D*Ts))/Ts + integrate(Y2,(t,0,(1-D)*Ts))/Ts

return average_values

In [34]: N(fun(0.5,*(400**2/75e3,100e-6,16e-6,200,10e-6 )),4)

Out[34]: 374.1
[ ]
399.2

In [35]: N(fun_r(0.5,*(400**2/75e3,100e-6,16e-6,200,10e-6 )),4)

Out[35]: 366.8
[ ]
391.4