Optimal Power Distribution Control for a Network of Fuel
Cell Stacks

arXiv:1601.07275v1 [] 27 Jan 2016

Resmi Suresh M P a, Ganesh Sankaran b, Sreeram Joopudi b,
Shankar Narasimhan a, Suman Roy Choudhury c, Raghunathan Rengaswamy a
a

Indian Institute of Technology, Madras, Chennai, Tamil Nadu, India 600036.
b
c

Gyan Data Private Limited, Taramani, Chennai, India

Naval Materials Research Laboratory, Ambernath, India

Abstract
In power networks where multiple fuel cell stacks are employed to deliver the required power, optimal sharing of the power
demand between different stacks is an important problem. This is because the total current collectively produced by all
the stacks is directly proportional to the fuel utilization, through stoichiometry. As a result, one would like to produce the
required power while minimizing the total current produced. In this paper, an optimization formulation is proposed for this
power distribution control problem. An algorithm that identifies the globally optimal solution for this problem is developed.
Through an analysis of the KKT conditions, the solution to the optimization problem is decomposed into off-line and on-line
computations. The on-line computations reduce to simple equation solving. For an application with a specific v-i function
derived from data, we show that analytical solutions exist for on-line computations. We also discuss the wider applicability of
the proposed approach for similar problems in other domains.
Key words: Power sharing; Optimization; Fuel cell; Analytical solution

1

Introduction

Fuel cells have gained considerable attention in the field
of power conversion, especially in the past few decades
[9,2,10]. Fuel cells are capable of extracting energy from
the fuel more efficiently than internal combustion engines [2]. Due to this higher fuel-efficiency, less harmful
emissions, low operational noise and absence of moving
parts, fuel cells have been proposed for a wide variety
of applications [7,2]. The power obtainable from a single
fuel cell is rather low in comparison to the power requirement for most of the applications. Hence, to meet realistic power requirements, active surface area has to be
increased. In view of this, multiple fuel cells are stacked
or connected in series and/or parallel arrangements to
enhance the output voltage and power [3,1].
While fuel cell technology has several advantages, there
are also problems related to long term durability and
gradual loss of performance over time. The performance
Email address: raghur@iitm.ac.in (Raghunathan
Rengaswamy).

Preprint submitted to Automatica

degradation could be due to many factors such as electrolyte degradation, catalyst poisoning, water flooding
etc. [2,3]. Such degradation can result in reduced current being drawn from the stack at the same voltage. As
a result, the overall system performance could be compromised, and in the worst case of complete stack failure, the whole system could fail [3]. This issue is usually
addressed by designing a reliable parallel connected fuel
cell architecture with an additional power bus as shown
in Fig. 1 [8]. Another approach to improve reliability is
the use of multiple power sources in the same network
[4,5]. While considerable research has been devoted to
concerns related to output voltage, state of charge, reliability and network configuration, very little work has
focused on developing operational strategies that minimize the fuel consumption through online control. We
address this control problem in this paper.
2

Problem Statement

Consider a fuel cell network with N branches and Nf,i
stacks in the ith branch (as shown in Fig. 2). The control
problem is one of minimizing the total current drawn (a

5 February 2018

branch is the decision variable. Subscripts lb and ub denote the lower and upper bounds respectively.
Assumptions:
1 The potential across each fuel cell is assumed to be a
static map fij :
Vij = fij (Ii )

∀ i = 1 : N;

j = 1 : Nf,i (2)

2 Number of stacks in each branch of the network under
consideration, is assumed to be one (Nf,i = 1 ∀ i =
1 : N ).
3 Minimum possible power from the fuel cell network
is achieved when all stacks are operating at their corresponding Ilb . Pathological cases where Pij |Ii,ub <
Pij |Ii,lb for j th stack in ith , which may arise depending
on the voltage-current profile, are neglected.

Fig. 1. Schematic of a parallel connected fuel cell network
with n branches

surrogate for total fuel utilized) while meeting a fixed
power requirement. This can be posed as the following
optimization problem.

Using equation (2), voltage constraints can be transformed to current constraints. Also, a network with multiple stacks in a branch can be converted to an equivalent network with a single stack in each branch, as current flow in each stack of a branch is equal. Power in
each branch of the resulting equivalent network can be
written as
Pi = Pef f,i = φi Vi Ii = φi fi (Ii )Ii

(3)

where Vi , Pi and φi are the voltage, power deliverable
and efficiency of the single fuel cell stack in the ith branch
of the equivalent circuit.
Using all these assumptions, the general constrained optimization problem described in (P1) can be reduced to
the following minimization problem (P2) for a network
as shown in Fig. 1.
Fig. 2. Schematic of a parallel connected fuel cell network
with N branches and Nf,i fuel cell stacks in the ith branch.

min
Ii

min
Ii

subject to

Inet =

N
X

Ii

subject to

(1a)

Pij = Preq

Vij ≤ Vij,ub

Ii

N
X

Pi =

N
X

(1b)

(4a)
φi fi (Ii )Ii = Preq

(4b)

i=1

Ii,lb − Ii ≤ 0
Ii − Ii,ub ≤ 0

i=1 j=1

Vij ≥ Vij,lb

N
X
i=1

i=1

i=1

f,i
N N
X
X

Inet =

∀ i=1:N
∀ i=1:N

(4c)
(4d)

The optimization problem (P2) in equation (4) involves
a linear objective function with one nonlinear equality
constraint and 2N linear inequality constraints. Due to
non-linear equality constraint present,this optimization
problem becomes non-convex.

∀ i = 1 : N ; j = 1 : Nf,i

(1c)
∀ i = 1 : N ; j = 1 : Nf,i
(1d)

where Inet is the total current from all branches, Ii represents the current drawn from the ith branch of the fuel
cell network. Vij , Pij (= φij Vij Iij ) and φij are voltage,
power produced and efficiency of the j th stack in the
ith branch respectively. Current (Ii ) drawn from each

3

Theoretical Analysis using KKT Conditions

Theorem 1 When a fuel cell network is used to achieve
a power of Preq , optimum is such that stacks operating at

2

Using

a current value within allowable limits (Ii,lb < Ii < Ii,ub )
i
are maintained at equal dP
dIi .

1−λ

Proof: For any V-I characteristics, Karush-KuhnTucker (KKT) conditions [6] can be used to prove the
above theorem. Lagrangian of the minimization problem described in (4) can be written as follows.

L=

N
X
i=1

Ii + λ Preq −

N
X
i=1

Pi

!

+

+

N
X
i=1

N
X
i=1

= 0,



dPi
dIi



γi (Ii − Ii,ub ) (5)

There are 2N inequality equations out of which first N
are lower bounds on Ii and the last N are upper bounds
on Ii . Necessary conditions for optimality as described
by KKT conditions for all i ǫ {1, ..., N } are

1−λ

(6)
(7)
(8)
(9)

1−λ
1−λ







dPi
dIi
dPi
dIi



For 1 ≤ M ≤ N and 1 ≤ K ≤ N ,
Constraints 1 : M are inactive, M + 1 : N are active.
Constraints N + 1 : N + K are inactive, N + K + 1 : 2N
are active.

=⇒

(a) If K ≤ M , for i = M + 1 : N , Ii = Ii,lb = Ii,ub
which is not feasible.
(b) If K > M , for i = K + 1 : N , Ii = Ii,lb = Ii,ub which
is not feasible.



∀

i = 1 + K : M (12)

− µi = 0 ∀ i = 1 + M : N (13)

Ii,lb



1


dPi
dIi

∀ i=K +1:M

(14)

Ii,opt



dPk+1
dIk+1



γj = −1 + 

(a) If K > M , For i = 1 + M : K, Ii = Ii,lb = Ii,ub
which is not feasible.
(b) If K ≤ M ,

for i = 1 + M : N

∀ i = 1 : K (11)

=

Ik+1,opt


dPk+2
= ...
dIk+2 Ik+2,opt


dPM
(15)
=
dIM IM,opt


i
Hence, dP
dIi of all stacks operating within their allowable limits are equal. Also, from (11) and (13),

For 1 ≤ M ≤ N and 1 ≤ K ≤ N ,
Constraints 1 : M are inactive, M + 1 : N are active.
Constraints N + 1 : N + K are active, N + K + 1 : 2N
are inactive.

Ii,lb

=0



Case 2

for i = 1 : K

(10)

Ii,opt is the optimum current drawn from ith
stack.From equation (12),

Case 1

Ii,ub

+ γi = 0

Ii,ub

Ii,opt

dPi
dIi

λ= 

Ii =

∀i=1:N

Substituting in equation (10),

Different cases in KKT conditions are obtained by varying the number of active and inactive constrains, all of
which falls into any of the following 2 cases for different
values of M and K.

(

− µi + γ i = 0

For the given case, KKT conditions can be reduced
to
(
0
for i = 1 : M
µi =
≥ 0 for i = 1 + M : N
(
≥0
for i = 1 : K
γi =
0
for i = 1 + K : N

µi (Ii,lb − Ii )

µi ≥ 0
γi ≥ 0
µi (Ii,lb − Ii ) = 0
γi (Ii − Ii,ub ) = 0

dL
dIi



µj = 1 − 

3

dPj
dIj
dPi
dIi

dPj
dIj

dPi
dIi









Ij,ub

∀ j = 1 : K;

Ii,opt

Ij,lb

i=K +1:M

(16)

∀ j = 1 + M : N;

Ii,opt

i=K +1:M

(17)

IN , then according to the algorithm, any incremental
power demand dPnet should be met by drawing additional current from those ’l’ such that,

It can be seen that γj ≥ 0 is satisfied if and only if


dPj
dIj



Ij,ub

≥



dPi
dIi



(18)

Ii,opt



and µj ≥ 0 can be achieved only if


dPj
dIj



Ij,lb

≤



dPi
dIi



dPj
dIj



Ij,lb

≤



≤





dPi
dIi

Ii,opt

Ii,opt

dPr
dIr

j = M + 1 : N; i = K + 1 : M;



(20)

Ir,ub

r = 1 : K (21)

Thus, first order optimality conditions imply that the
only feasible solution to the concerned optimization
problem (P2) is such that ’M’ stacks are operating at
a current higher than their minimum allowable current
(Ii,lb ) and out of which ’K’ stacks are operating at Ii,ub
and the conditions (15) and (20) are satisfied.
4

A new algorithm is proposed to find an optimum solution
to the minimization problem (P2) described using equation (4). For a network of ’N’ branches with a single stack
in each branch, to deliver Preq , say Ii > Ii,lb for first ’m’
stacks and among these, ’k’ stacks operate at Ii = Ii,ub
at optimum. Stacks k + 1 to m are then operating at a
current in between their bounds, Ii,lb < Ii < Ii,ub . At
optimum, according to theorem 1,
dPi
dIi



=

Ii,opt



dPj
dIj



5

subject to
m
X

ef f
Pi = Preq

(23a)

i=k+1



dPj
dIj



Ij,lb

≤



dPi
dIi



Ii,opt

such that j = m + 1 : N ;
and
ef f
Preq
= Preq −

N
X

j=m+1

≤



dPr
dIr



Pj,min −

k
X

(23b)

Ir,ub

i = k + 1 : m;

Pr,max

≥

dPnet
dIn



In

∀ nǫN

(25)

Implementation

The implementation of proposed algorithm is described
in Fig. 3. Input to the algorithm will be the V-I model
parameters, bounds on current and the required power
(Preq ). Using the new algorithm, majority of the calculations for estimating optimum current can be performed
off-line, reducing the on-line computations. The major
step in the off-line mode is to identify observable points
(kink points which decide whether to add or remove
stacks from the list of stacks, for which current needs to
be optimized). Observable points are decided based on
(dPj/dIj ) of each stack evaluated at both Ij,lb and Ij,ub ,
where j corresponds to the branch number. Each value in
the list of (dPj/dIj ) corresponds to an observable point.
For a fuel cell network with N branches, there will be 2N
observable points. These observable points are then arranged such that the power corresponding to the observable point decreases down the list. Corresponding to the
ith value of (dP/dI ) in the sorted list, if for any j ǫ 1 : N ,
(dPj/dIj )Ij =Ij,lb ≤ (dP/dI )i , then j th stack should be operating at Ij,lb . And for any j ǫ 1 : N , if (dPj/dIj )Ij =Ij,lb >
(dP/dI )i > (dPj/dIj )Ij =Ij,ub , then j th stack should be operating at a current in between Ij,lb and Ij,ub . Optimum current to be drawn from these stacks can be obtained by solving (dPj/dIj )Ij,opt = (dP/dI )i . Further, if
for any j ǫ 1 : N , (dPj/dIj )Ij =Ij,ub ≥ (dP/dI )i , then j th

∀ i, j = k + 1 : m (22)

Ij,opt

Il



Remark 1 The optimization problem described in (4) is
convex for any V-I profile, Vi = fi (Ii ), with a restriction
that corresponding power (Pi = φi Ii fi (Ii )) is a concave
function (or in other words Hessian of power is negative
definite). For this convex subproblem, though multiple
solutions are obtained by solving equations (22), (23) and
(24), only one of them will satisfy the actual constraints
of the optimization problem (P2), resulting in a unique
solution. Proof can be found in Appendix A.

Proposed Algorithm





For a power requirement of Preq , any configuration that
satisfies the above criterion will be a local optimum to
the problem. Multiple solutions are possible depending
on the V-I profile. Global optimum can be obtained
based on the value of objective function for those solutions which satisfy all these constraints. The solution
which has lowest value of objective function will be the
global optimum. This optimum will have maximum rise
in total power from
for a specific increase in
 thenetwork
PN dPi 
dPnet
= i=1 dIi .
current drawn
dInet

(19)

Combining (18) and (19),


dPnet
dIl

r=1:k

(24)

r=1

Corresponding to a power requirement of Preq , let the
optimal current drawn from branches be I1 , I2 ,..., and

4

on-line such that equations (22), (23) and (24) are satisfied. The exact solution to this set of equations has been
obtained analytically for a specific function Vi = fi (Ii )
and is discussed in section 6. If multiple solutions are
obtained by solving these set of equations which satisfy
the bounds on current, then global optimum is found
based on the value of objective function. The solution
with lowest value of net current is the global optimum
for the problem.

Start
f (Ii ), Iub , Ilb

Calculate Pobt,min and Pobt,max

Arrange cells in decreasing order of ( dP
dI )
calculated for each cell at Ii,lb and Ii,ub .
i = 1 : 2N represents index of sorted list
Compute Pobt at each observable points

6
Preq

Pobt,min ≤
Preq ≤ Pobt,max

no

Error: Required
power cannot
be obtained

Analytical solution for a specific convex subproblem
V−I characteristics of a fuel cell
26

Experimental Data
Model Fit

n=0 yes
n=n+1

no

Voltage (V)

Preq ≤ Pobt,n

24

Update n

yes Corresponding observable point = n-1
Find cells operating at Iub (k cells) and Ilb (N − m cells)

Model Obtained
V = 27.3873 − 0.7804 I0.5

22
20
18

Find cell numbers for which current
needs to be optimized (m − k cells)

16
0

Calculate
ef f
Preq
(24)

50

100
Current (A)

150

200

Fig. 4. Comparison of experimental data obtained from
NMRL with the assumed V-I characteristics form.

i
Solve ( dP
dIi ) = constant ∀ i = k + 1 : m
Pk
ef f
and
= Preq
i=1 Pi

Voltage and current data (See Fig. 4) for a typical fuel
cell stack are obtained from Naval Materials Research
Laboratory (NMRL). V-I characteristics of the form as
described in (26) is assumed for which corresponding
power is a concave function. The assumed profile has
been tested on the data and a 99.87% good fit is obtained
as shown in Fig. 4.

Multiple solutions

Current in cells i = k + 1 : m
Current ∀ i = 1 : k & i = m + 1 : N
All possible local optima
Minimum total current
Optimum
Current: Iopt
Stop

p
V (Ii ) = fi (Ii ) = ai + bi Ii
such that ai ≥ 0
bi ≤ 0
∀ i=1:N

Online mode

Fig. 3. Proposed algorithm: A fuel cell network with N
branches is considered such that, to deliver the required
power (Pr eq) m cells should be operated at a current higher
than Ilb , out of which k cells operate at Iub .

(26a)
(26b)
(26c)

As power corresponding to this V-I profile is concave,
corresponding optimization problem is convex as described in Appendix A and thus will have a unique global
solution. Analytical solution is derived for this specific
V-I characteristics.

stack should be operating at Ij,ub . Hence, corresponding to each observable point, current in each stack is
pre-calculated off-line. This information can be used to
evaluate the power at every observable point.

Consider ’n1’ cells operating at Ij,opt at optimum such
that Ij,lb < Ij,opt < Ij,ub . The following equations can
be obtained using Theorem 1 and are used to determine
optimum current drawn from those ’n1’ cells as shown
in Fig. 3.
ef f
P1 + P2 + ... + Pn1 = Preq
(27)

When there is a specific power requirement, the stacks
to be operated at current values between their bounds
are identified based on the pre-calculated power at each
observable point. Upper and lower bound constraints on
current are inactive for these stacks. The optimum current to be drawn from each of these stacks is evaluated

dP1
dP2
dPn1
=
= ... =
dI1
dI2
dIn1

5

(28)

Using equation (3) and (26), and substituting xj =
equation (27) and (28) reduces to
n1
X

ef f
(φj aj + φj bj xj ) x2j = Preq

p
Ij ,

multiple stacks in a branch. Finally, efficiency of the
algorithm is compared with that of standard optimizers.
7.1

(29)

j=1

 


3
3
φ1 a1 + φ1 b1 x1 = φj aj + φj bj xj
2
2

∀ j = 2 : n1
(30)

For j = 1 : n1, rearranging equation (30),
xj = gj x1 + hj
= gj2 x21 + 2gj hj x1 + h2j


3
3 3
xj = gj x1 + x21 3gj2 hj + x1 3gj h2j + h3j
x2j

(31)
(32)
(33)

(φ1 a1 − φj aj )
and hj =
1.5φj bj

(34)

Substituting in equation (29) and simplifying, we get a
cubic equation.

x31




n1
X
j=1



φj bj gj3 

+

x21




n1
X
j=1

(35)

j=1

Current drawn from all stacks can be evaluated using
equations (31) and (36) after solving for x1 (using equation (35)).
Iopt,j = x2j
(36)
As the analytical solution is obtained by solving a cubic
equation (Equation (35)), it will result in three sets of
solutions. But as the problem is convex and has a unique
minimum, it is guaranteed that there can only be one
solution set among those three which will satisfy all constraints of the actual optimization problem (P2).
7

Stack

Lower bound

Upper bound

Number

on current (A)

on current (A)

φ

a

b

1

2.103

2

0

106.8127

1

47.655

-1.297

325.6562

1

39.895

-0.557

3

6.646

236.4155

1

33.847

-0.5976

stacks are operating at their corresponding lower bounds
is 310.976 W. For each observable point, stacks which
are to be operated at their corresponding lower bounds
and upper bounds can be pre-calculated and are given in
Table 2. It can be seen from both Table 2 and Fig. 5 that
stack 1 starts operating first on increasing the power requirement from Pmin (= 310.976W ). When the power
required reaches the power corresponding to the second
observable point (corresponding to dP2/dI2 evaluated at
I2,lb ), then stack 2 also starts to operate. Each stack
starts to operate one by one and for the stacks operating, optimum currents are calculated using analytical
solution described in section 6. If observable point corresponding to upper bound current of stack j is reached,
then from that power onwards ’j th ’ stack should be operating at its upper bound current. This can be seen starting from observable point 4, where stacks reaches their
corresponding upper bounds one by one. The process
continues in a similar manner till all stacks are operating
at their upper bound. Fig. 5 shows how various stacks
should be operated at optima for various power requirement regimes. Fig. 6 shows the optimum current calculated using the proposed algorithm for various power
requirements.



3φj bj gj2 hj + φj aj gj2 



n1
X

ef f
3φj bj gj h2j + 2φj aj gj hj  = Preq
+ x1 

Consider an example of a fuel cell network consisting of
3 fuel cell stacks that are connected in parallel using a
power bus. Parameters used for simulation√are described
in Table 1. The variation of dPi/dIi with Ii for all the
stacks are shown in Fig. 5. It can be clearly seen that
stack 1 has the highest dPi/dIi at its lower bound and so it
should start operating first according to the algorithm.
The observable points as explained in section 4 are calculated and tabulated in Table 2 and are also marked
in Fig. 5. Power obtainable from the network when all
Table 1
Parameters for a fuel cell network with 3 fuel cell stacks

where

φ1 b1
gj =
φj bj

Working of algorithm

Results and Discussion

The algorithm proposed in section 4 has been implemented on various fuel cell networks. All fuel cell stacks
are assumed to follow V-I characteristics described in
equation (26). Parameters for various stacks are obtained by perturbing the parameters calculated for a
typical fuel cell stack V-I profile shown in Fig. 4. The
results obtained for three different case studies are presented here. The working of algorithm is shown using
the first case study while the second one aims at demonstrating the application of algorithm for a network with

For any power requirement, we get three solutions by
solving equation (35). For example, when Preq = 8000,
three solutions obtained are given in Table 3. But as explained in section 6, it can be seen that only one solution
satisfies upper and lower bounds for current. It is an interesting inference that only a single solution will satisfy
all the constraints out of all the solutions obtained using
algorithm for a convex problem.

6

Fig. 5. dPi/dIi variation
for different stacks in the network for which parameters are described in Table 1. Variation is shown
√
with respect to I such that I varies from Ilb to Iub . Each horizontal line represents an observable point.
Table 2
Observable points calculated for the fuel cell network with 4 stacks whose parameters are described in Table 1
Sl.no

Corresponding

Corresponding

Stacks operating

Stacks

Stacks operating

total power

dP/dI

at Imin

to be optimized

at Imax

1

310.976

44.834

1,2,3

-

-

2

890.577

39.895

2,3

1

-

3

6183.777

31.536

3

1,2

-

4

12037.033

27.548

-

2,3

1

5

16200.563

24.818

-

3

1,2

6

19206.708

20.064

-

-

1,2,3

Table 3
Three solutions obtained by solving equation for the fuel
cell network with parameters shown in Table 1 for a power
requirement of 8000 W.
Stack 1

Stack 2

Stack 3

Total Current

Solution 1

1140.37 A

4808.83 A

3350.97 A

9300.17 A

Solution 2

81.02 A

136.23 A

17.07 A

234.32 A

Solution 3

1.93 A

36.6 A

153.407 A

191.94 A

7.2

ble 4. For applying the algorithm, an equivalent network
with single stack in each branch is constructed for this
network such that total power from each branch of the
original network is the same as the power from that single stack in the corresponding branch of the equivalent
network.

Let Pi be the power acquired from ith branch in the
equivalent network, then

Equivalent network formation for a network with
multiple stacks in each branch

Nf,i

Nf,i

Consider a network of 30 fuel cell stacks unequally distributed in 15 branches with parameters as given in Ta-

Pi =

X
j=1

7

Pij =

X
j=1

φij Vij Ii

(37)

Table 5
Comparision of the results obtained for unequal distribution
of fuel cell stacks in different rows using a standard optimization routine and the proposed algorithm. Preq =75000
W, φ=0.8 Number of fuel cell stacks = 30. Parameters and
optimum current drawn are described in Table 4.

350
Cell 1
Cell 2
Cell 3

300

Current (A)

250

Optimizer

Proposed Algorithm

Total Current (A)

828.88

828.88

Total Power (W)

75000

75000

0.104

0.0024

200
150
100

Time taken for optimizing (sec)

50

sets of solutions obtained from the algorithm satisfied
upper and lower bound constraints. It is surprising to
see that a cubic equation as given in (35) will always result in one and only one solution that satisfies all the
constraints for a network with all stacks having a concave power profile. For the convex subproblem, a unique
optimum is guaranteed from the algorithm.

0
0

0.5

1
Power (W)

1.5

2
4

×10

Fig. 6. Optimum currents obtained for stacks for various
power requirement

For a V-I profile given in equation (26),

For a general problem, if required power is feasible from
the fuel cell network, then all possible local minima can
be obtained from the proposed algorithm unlike conventional optimizers where a single local optimum is obtained based on the initial guess. Also, all these possible local optima are estimated in a very short time.
Global optimum is guaranteed using the proposed algorithm unlike standard optimizers as it can be estimated
easily using all possible local optima. If power is a concave function of current for all cells in the network, then
a unique local optimum is obtained from the algorithm
and this local optimum will be the global optimum to
the problem. Optimum solution and computational time
are highly dependent on the initial guess fed by the operator in case of any optimization routine. Whereas, the
algorithm proposed does not require an initial guess and
is much more computationally efficient.

Nf,i

Pi =

X

p
φij (aij + bij Ii )Ii

(38)

X

(φij aij Ii + φij bij Ii1.5 )

(39)

j=1

Nf,i

=

j=1

= (ai + bi

p
Ii )Ii

(40)

PNf,i
PNf,i
where ai =
j=1 (φij aij ) and bi =
j=1 (φij bij ) are
parameters in voltage profile for the stack in ith branch of
the equivalent network. Using the equivalent network of
the above mentioned network, optimum current in each
branch is estimated for a power requirement of 75000 W
and is given in Table 4.
7.3

Comparison of proposed algorithm with conventional optimizers

Though this problem can be solved using standard optimizers in MATLAB or Python, it is not possible to embed this optimizer in a chip due to various limitations
given below.

The optimum current estimated in each branch of various fuel cell networks using proposed algorithm is compared with that of standard optimizers. The optimum
currents obtained using both approaches are found to
be matching well for all the networks. Table 4 shows the
comparison of optimum current estimated using standard optimizer and proposed algorithm in Python environment for one of the networks. The time taken for
optimization and the value of objective function using
both approaches when implemented in Python are given
in Table 5. It can be seen that the proposed algorithm
is much faster than the standard optimizer.

(a) It cannot determine whether a solution to the problem
is feasible or not.
(b) It is sensitive to the initial guess passed to the optimizer.
(c) Robustness cannot be ensured.
(d) Increased computational load on the hardware.
(e) It requires large computational time and storage
space.
(f) It does not guarantee global solution.

Multiple trials were run for different power requirements
and for various fuel cell networks where all stacks are assumed to have a concave power profile.. It was observed
that in all these runs, only one solution among the three

All the above limitations can be overruled by the new
algorithm and thus can be embedded into a chip for online applications.

8

Table 4
Optimum current values obtained in case of unequal distribution of fuel cell stacks in different rows. Preq =75000 W, φ=0.8,
Number of Fuel Cell Stacks = 30, 0.1 ≤ I ≤ inf

Stack

Row

a

b

Iopt (A)

Iopt (A)

No.

No.

values

values

(Optimization routine)

(Proposed Algorithm)

1

1

49.25

-0.25

2

1

49.302

-0.302

9.5299

9.5355

3

2

49.353

-0.353

0.1

0.1

4

3

49.405

-0.405

5

3

49.457

-0.457

648.8999

648.9024

6

3

49.509

-0.509

7

4

49.56

-0.56

8

4

49.612

-0.612

3.2599

3.262

9

5

49.664

-0.664

10

5

49.716

-0.716

2.6699

2.674

11

6

49.767

-0.767

12

6

49.819

-0.819

2.2799

2.2766

13

7

49.871

-0.871

14

7

49.922

-0.922

1.99

1.993

15

8

49.974

-0.974

0.1

0.1

16

9

50.026

-1.026

17

9

50.078

-1.078

1.6899

1.6946

18

10

50.129

-1.129

19

10

50.181

-1.181

1.55

1.55

20

11

50.233

-1.233

21

11

50.284

-1.284

90.0699

90.0715

22

11

50.336

-1.336

23

12

50.388

-1.388

24

12

50.44

-1.44

1.3

1.2976

25

13

50.491

-1.491

26

13

50.543

-1.543

64.1899

64.1877

27

13

50.595

-1.595

28

14

50.647

-1.647

29

14

50.698

-1.698

1.13999

1.1365

30

15

50.75

-1.75

0.1

0.1

9

8

Conclusions

[7] James Larminie, Andrew Dicks, and Maurice S McDonald.
Fuel cell systems explained, volume 2. Wiley New York, 2003.
[8] Steven C Michalske and Bradley L Spare. Parallel fuel
cell stack architecture. April 15 2010. US Patent App.
12/761,284.

An optimization problem is formulated for the optimum
power distribution control for a fuel cell network. A new
and computationally efficient algorithm is proposed for
this optimization problem which can easily be embedded
in a chip for online applications. New algorithm is proposed based on the observations when KKT conditions
are applied to the given optimization problem. The algorithm proposed here guarantees global optimum to the
optimization problem concerned unlike local optimum
obtained using the conventional method of using standard optimizers and in a much shorter time. This can be
very useful in optimizing the current and implementing
control strategies in a real fuel cell network. Analytical
solution for a specific convex subproblem is also formulated in this work. Though the analytical solution is applicable only for a specific voltage current profile, the algorithm can be used for any other V-I characteristics as
well, irrespective of whether the resulting optimization
problem is convex or not.

[9] B Rohland, J Nitsch, and H Wendt. Hydrogen and fuel
cellsthe clean energy system. Journal of power sources,
37(1):271–277, 1992.
[10] Colleen Spiegel. PEM fuel cell modeling and simulation using
MATLAB. Academic Press, 2011.

A

Consider the optimization problem (P3) where the
equality constraint is relaxed to inequality constraint.
min

Inet =

Ii

N
X

Ii

(A.1a)

i=1

subject to

The proposed algorithm can be further extended to other
power distribution systems including pumps and compressors, generators, batteries and solar cells, though a
fuel cell network is considered through out this paper.
For example, in a multiple pump system used for water
transport with varying demand, the algorithm can determine which pump should operate at what flow rate to
satisfy the required water demand at any point of time.
Here, flow rate is analogous to current while head or
pressure drop is comparable to voltage. Also, electrical
systems used in a parallel network configuration is similar to the fuel cell network and its optimum current distribution among connected sources, to achieve a desired
power, can be identified using the algorithm. It can also
be applicable for cases where combinations of different
systems are connected together to form a network.

N
X
i=1

Pi ≥ Preq

(A.1b)

Ii,lb − Ii ≤ 0

∀ i=1:N

(A.1c)

Ii − Ii,ub ≤ 0

∀ i=1:N

(A.1d)

Lagrangian for this problem is

L=

N
X
i=1

Ii + λ Preq −

N
X

Pi

i=1

+

!

+

N
X
i=1

N
X
i=1

µi (Ii,lb − Ii )

γi (Ii − Ii,ub )

(A.2)

There are a total of 2N+1 inequality constraints out
of which first N are lower bounds on current, next N
are upper bounds on current and the last one is the relaxed power inequalty constraint. Necessary conditions
for optimality as described by KKT conditions for all
i ǫ {1, ..., N } are

References
[1] How fuel cells work: Polymer exchange membrane fuel cells.
Accessed on 08 January 2015.
[2] L Carrette, KA Friedrich, and U1 Stimming. Fuel cells–
fundamentals and applications. Fuel cells, 1(1):5–39, 2001.
[3] Fuel Cell Handbook.
Eg&g technical services.
Albuquerque, NM, DOE/NETL-2004/1206, 2004.

Relaxed Optimization Problem

µi ≥ 0
γi ≥ 0
µi (Ii,lb − Ii ) = 0
γi (Ii − Ii,ub ) = 0
λ≥0

Inc.,

[4] Zhenhua Jiang and Roger A Dougal. Strategy for active
power sharing in a fuel-cell-powered charging station for
advanced technology batteries.
In Power Electronics
Specialist Conference, 2003. PESC’03. 2003 IEEE 34th
Annual, volume 1, pages 81–87. IEEE, 2003.

λ(Preq −

[5] Zhenhua Jiang, Lijun Gao, and Roger A Dougal. Adaptive
control strategy for active power sharing in hybrid fuel
cell/battery power sources.
Energy conversion, ieee
transactions on, 22(2):507–515, 2007.

N
X

Pi ) = 0

(A.3)
(A.4)
(A.5)
(A.6)
(A.7)
(A.8)

i=1

We get 2 additional conditions using the relaxed constraint. Different cases using KKT conditions can be
classified to mainly 2 cases as given below.

[6] HW Kuhn and AW Tucker. Proceedings of 2nd berkeley
symposium, 1951.

10

as the cases given in section 3 and we get the same conditions for feasible solution as the original optimization
problem (P2) given by

Case I
Power constraint is inactive: Preq <

PN

i=1

Pi and λ = 0.



1) For 1 ≤ M ≤ N and 1 ≤ K ≤ N ,
Constraints 1 : M are inactive, M + 1 : N are active.
Constraints N + 1 : N + K are inactive, N + K + 1 : 2N
are active.

dPk+1
dIk+1



=

Ik+1,opt



dPk+2
dIk+2



= ...

Ik+2,opt

=





dPj
dIj



Ij,lb

≤



dPi
dIi



Ii,opt

≤



j = M + 1 : N; i = K + 1 : M;
2) For 1 ≤ M ≤ N and 1 ≤ K ≤ N ,
Constraints 1 : M are inactive, M + 1 : N are active.
Constraints N + 1 : N + K are active, N + K + 1 : 2N
are inactive.

Using

dL
dIi

Ii,ub

for i = 1 : K

Ii,lb

for i = 1 + M : N

A.1

∀i=1:N

(A.9)

Substituting in equation (A.9),
∀ i=1:K
∀ i=1+K :M
∀ i=1+M :N



(A.14)

Ir,ub

r=1:K

(A.15)

Convex Subproblem

The solution to this convex problem can be obtained
from the same conditions derived in section 3. But, as it
is a convex problem, only one minima is possible while
equations 22 and 23 will result in multiple solutions.
Hence, among all possible multiple solutions to those
conditions, only one of the solutions will satisfy the
bounds on current and results in a single optima to the
convex subproblem.

For the given case, KKT conditions can be reduced
to
(
0
for i = 1 : M
µi =
≥ 0 for i = 1 + M : N
(
≥0
for i = 1 : K
γi =
0
for i = 1 + K : N

1 + γi = 0
1=0
1 − µi = 0

dPr
dIr

Let the power of each stack (Pi = φi Ii f (Ii )) be a concave function (or in other words Hessian of power be
negative definite), then all the constraints in optimization problem (P3) given by equation (A.1) together form
a convex set. Moreover, the objective function is linear.
Thus, the minimization problem (P3) is a convex problem for a concave power function.

= 0,

1 − µi + γ i = 0

(A.13)

IM,opt

Hence, the relaxed optimization problem has a feasible
solution if and only if the power constraint is active.
Moreover, the solution set is the same for both minimization problems P2 and P3. Thus, those two optimization
problems can be considered equivalent.

(a) If K > M , For i = 1 + M : K, Ii = Ii,lb = Ii,ub
which is not feasible.
(b) If K ≤ M ,
Ii =



and

(a) If K ≤ M , for i = M + 1 : N , Ii = Ii,lb = Ii,ub
which is not feasible.
(b) If K > M , for i = K + 1 : N , Ii = Ii,lb = Ii,ub which
is not feasible.

(

dPM
dIM

(A.10)
(A.11)
(A.12)

which cannot be satisfied. Hence, this case is also not a
feasible solution.
Case II
PN
Power constraint is active: Preq = i=1 Pi and λ ≥ 0.
The sub cases when power constraint is active is same
11

Observable points

Current I (A)

20

400

Obs. pt. 1
Obs. pt. 2
Obs. pt. 3
Obs. pt. 4

350
Current I (A)

25

15

10

5

0
1

Obs. pt. 5
Obs. pt. 6
Obs. pt. 7
Obs. pt. 8

300

250

200

2
3
Cell Number

4

150
1

2
3
Cell Number

4

Variation of dP/dI with I0.5
45
Cell 1
Cell 2
Cell 3
Cell 4

(W/A)

35

dP/dI

40

30

25

20
0

5
0.5

I

10
(Ampere0.5)

15

20

Optimum currents in each cell for various power requirement
400

P =100
P =350
P =1000
P =15000
P =31160
P =33000
P =34800

350

Current I (A)

300
250
200
150
100
50
0
1

1.5

2

2.5
Cell Number

3

3.5

4

