Post-Prognostics Decision for Optimizing the Commitment of Fuel
Cell Systems
Stéphane Chrétiena,1,∗, Nathalie Herrb , Jean-Marc Nicodb , Christophe Varnierb
a

arXiv:1710.08812v1 [] 19 Oct 2017

b

National Physical Laboratory
FEMTO-ST/AS2M, Univ. Bourgogne Franche-Comté/UFC/CNRS/ENSMM,Besan con, France

Abstract
In a post-prognostics decision context, this paper addresses the problem of maximizing the
useful life of a platform composed of several parallel machines under service constraint. Application on multi-stack fuel cell systems is considered. In order to propose a solution to
the insufficient durability of fuel cells, the purpose is to define a commitment strategy by
determining at each time the contribution of each fuel cell stack to the global output so as
to satisfy the demand as long as possible. A relaxed version of the problem is introduced,
which makes it potentially solvable for very large instances. Results based on computational
experiments illustrate the efficiency of the new approach, based on the Mirror Prox algorithm, when compared with a simple method of successive projections onto the constraint
sets associated with the problem.
Keywords: Decision making, Post-prognostics decision, PHM, Fuel cell, Convex
optimization example, LATEX, template

1. Introduction and related work
In the context of the decline of fossil fuel resources, the use of fuel cells appears to be
of growing interest as a potential alternative to conventional power systems [10]. Fuel cells
can be used in many applications, such as stationary ones for domestic use, but also in
transportation and portable power applications [4]. Unfortunately, fuel cells suffer from
insufficient durability. Indeed, the current lifetime is usually between 1 500 and 3 000 hours,
whereas 5 000 hours are often required for transportation applications and up to 100 000
hours for stationary ones. Improving the performance, reliability and lifetime is therefore
an important challenge [4] in the fuel cell technology, for which techniques of Prognostics
and Health Management (PHM ) can definitely help. As recently pointed out in Jouin et al.
∗

Corresponding author
Email addresses: stephane.chretien@npl.co.uk (Stéphane Chrétien), nathalie.herr@femto-st.fr
(Nathalie Herr), jean-marc.nicod@femto-st.fr (Jean-Marc Nicod), christophe.varnier@femto-st.fr
(Christophe Varnier)
URL: https://sites.google.com/view/stephanechretien/home (Stéphane Chrétien)
Preprint submitted to Computers & Operations research

October 25, 2017

[10], research in fuel cells with a PHM viewpoint has mainly focused on data acquisition
and data processing. Less work has been done on condition assessment and diagnostics
and very few papers addressed prognostics and decision making. Works that take into
account the decision part often only proposed corrective actions (see [5] and [14]), for which
physical parameters (such as inlet and outlet gas flows, pressures and temperatures, single
cell and stacks voltages or current) are controlled in order to master each fuel cell operating
conditions as accurately as possible. These corrective actions correspond to real-time control
(from nanoseconds to seconds), necessary to compensate the natural fluctuation of fuel cells
parameters and to avoid too early irreversible degradation. At each time slot the strategy
allows to set the operating current to a value that meets the needs in power for each fuel
cell.
Decision making as addressed in this paper differs from the previous studies in many
aspects. In particular, larger time scales (hours to weeks) are considered and decision comes
within the framework of Prognostic Decision Making (PDM), which aims at optimizing
systems configuration [1]. The system considered here consists of several fuel cells, used in
parallel in order to provide a global power output. The problem we consider is to provide
the power output value for each fuel cell as a function of time, so as to meet the global power
demand. The target application considered here is based on stationary power generation for
domestic usage, also known as micro combined heat and power (micro-CHP).
In order to deliver suitable power outputs, fuel cells are used in the form of stacks,
composed of many connected single cells. Each stack is assumed to be independent of the
others, but the multi-stack fuel cell system has to deliver a given global power output. At
each time slot, the total provided power output is the sum of each output of the stacks that
are switched on. Each fuel cell stack is able to deliver an output that can vary continuously
within a given interval. The optimization problem consists in determining the appropriate
output for each fuel cell stack during the whole production horizon. Also the stacks are
not constrained to be all running at each time slot if the target output can be reached by
using only a subset of them. All the stacks may moreover not be available at all times:
if their end of life has been reached, they are not available for production. Considering
a global power output, the multi-stack system useful life depends not only on each stack
useful life, but also on both the schedule and the operating condition settings that define
the contribution of each stack over time. The same statement applies to batteries in a
health management context. Saha et al. [13] have for instance addressed the maximization
of the battery charge used while constraining the probability of a battery shut off in flight for
electric unmanned aerial vehicles. Predictions on remaining battery life are used to optimize
mission plans without exceeding the available battery charge. In a same way, we propose to
use prognostics results in the form of Remaining Useful Lives (RUL) to maximize the global
useful life of a multi-stack fuel cell system subject to service constraint.
A similar problem has been addressed in [7] and [8], where the purpose was to define
a schedule of machines that maximizes the production horizon, based on the knowledge of
each machine remaining useful life (RUL) in a PHM framework. In these studies, machine
throughputs have been assumed to take discrete values. It was shown in [8] that optimal
solutions can be found in reasonable time only for small size instances and with a very limited
2

number of machines, very few possible throughput values and short production horizons.
An other study considering this time machines whose performances can vary continuously
between two bounds has been proposed in [9]. The considered model has been built to fit
the fuel cells behavior, but the proposed resolution approach gives suboptimal solutions and
is limited to systems of reasonable size.
In order to overcome these two limitations, we propose in this paper to drastically change
the solution paradigm and to build the solutions globally on the full production horizon.
Contribution of each machine during its lifetime is considered as a whole and optimized
on the full horizon. Each machine contribution is determined via convex optimization. As
a result, our new approach allows to address large problems very quickly. Concretely, the
considered scheduling problem is addressed via optimizing a composite function subject to
several constraints due to fuel cell intrinsic characteristics. A first method is described
which performs successive projections onto the sets of constraints. The method generates
a sequence of points that can be shown to converge to the solution of the optimization
problem [2]. A second order method is developed and used to define the contribution of
each fuel cell stack to the global output over the whole production horizon. It is based on
the Mirror Prox method proposed by Nemirovski [11] as a variant of the Mirror Descent
developed by Nemirovsky and Yudin [12] to minimize a smooth convex function subject to
convex constraints. Estimation of the variable is efficient in that it depends very little on its
dimension. This is why these methods can be used to solve big optimization problems [3].
The organization of the paper is as follows: the tackled problem is first described in
Section 2, with a presentation of the application framework and the optimization problem,
followed by a mathematical formulation that complies with the proposed convex resolution
paradigm. The two proposed resolution approaches are then described in Sections 3 and 4.
Efficiency of these methods is assessed through simulation results in Section 5. Conclusion
and future work are finally given in Section 6.
2. Problem statement
2.1. Application framework
The application addressed in this paper is based on a multi-stack fuel cell system which
is supposed to meet energy requirements for domestic usage in a stationary power generation
framework. This system is supposed to be composed of m fuel cell stacks Mj (1 6 j 6 m).
All the stacks are supposed to be always supplied with raw material required for the energy
conversion. They can be used simultaneously and independently from each other.
This corresponds to a parallel machines system, in which each machine is supposed to
be able to deliver power outputs Pj that can vary continuously within a given power output
range [P minj ; P maxj ]. For each machine Mj (1 6 j 6 m), the range of available power
outputs depends on the time t and evolves as depicted in Figure 1, with the following extreme
output characteristics:
P maxj : the maximal power output which decreases with time when the machine Mj is used.
The straight line depicting its decrease has equation P maxj (t) = aj · t+P maxj (0), with
3

aj < 0 the speed of the output decline and P maxj (0) the maximal output available
at the beginning of the scheduling process. Both parameters of this equation are fuel
cells intrinsic characteristics. Useful life of each power output Pj , denoted RULj (Pj ),
is limited by this maximal power output decrease ;
P minj : the minimal power output, constant over time and associated to the maximal Remaining Useful Life: RULj (P minj ) = maxP minj 6Pj 6P maxj (RULj (Pj )) = RULmaxj .
Pj (t)
P maxj (0)

P maxj (0)
100%

P maxj (t) = aj · t + P maxj (0)
Pj

Decrease
of P maxj
P minj
0

15%

RULmaxj

RULj (Pj )

time

(End of life)

Figure 1: Evolution of available power outputs for a machine Mj

The model proposed in Figure 1 is a simplified version of the one developed in [9]. It
considers however the main characteristics of fuel cells useful for the scheduling purpose
considered here. Due to the decrease over time of the maximal power output P maxj (t), this
model is asymmetric. In a schedule defining the contribution over time of a machine, two
periods of time in which a machine is used with different outputs can then not always be
permuted (see [9] for an illustration of this property).
2.2. Optimization problem
At each time t (0 6 t 6 T ), the global outcome P tot is the sum of each stack power
contribution. During the whole production horizon, this global outcome has to reach a
given load demand σ(t). Overproduction is authorized if it allows to extend the global
system useful life (or production horizon), but should be avoided as far as possible. In
fact, storage being not considered in this study, overproduction is supposed to be lost. All
the machines are not supposed to be in use at each time if a subset of them is enough to
reach the demand. Since machines suffer from wear and tear, some machines can also be
unavailable at a certain time if their end of life (EOL) has been reached. Stop-and-start of
fuel cell stacks have moreover to be avoided as far as possible. Stopping and restarting a
fuel cell can indeed induce considerable damage and lead to premature aging [4]. Change of
power output during the use of fuel cell stacks is however still authorized.
Considering these assumptions, the point is to manage the system by defining the commitment of fuel cell stacks so as to reach the demand as long as possible. During the whole
production horizon, the purpose is then to define at each time each stack contribution to
the global power output.
4

2.3. Mathematical formulation
In this section, a mathematical formulation of the problem previously described is defined, making use of convex elements. Let fj (t) (1 6 j 6 m, 0 6 t 6 T ) be the vector
defining the evolution over time of the power output delivered by the machine Mj , with
T the length of the decision horizon. Link between this decision horizon and the solution
production horizon, denoted H, is clarified in Section 5.3. Contributions of all the machines
are gathered together in a vector F ∈ Rm(T +1) such that:
F = [f1 (0), . . . , f1 (T ), . . . , fj (t), . . . fm (0), . . . , fm (T )] .
The general idea is to minimize a convex function subject to a set of constraints. The
objective function aims at ensuring that the power demand is reached. At each time t
(0 6 t 6 T ), it is about minimizing the difference between the global power output delivered
by the set of machines and the demand σ(t). This is expressed by Equation (1), where φ
measures the error incurred by the choice of fj , ∀ 1 6 j 6 m. One possible choice of φ is
the squared loss, i.e. the squared L2 norm.
!
m
X
min φ σ(·) −
fj (·)
(1)
j=1

Constraints on each function fj relate to the definition domain of each contribution and
to the limited availability of machines. At each time t, each machine contribution is either
equal to zero or constrained between two bounds (see Equation (2)), in accordance with the
hypotheses detailed in the application framework. fj (t) = 0 means that the machine Mj is
not used at time t.
fj (t) = 0 or fminj (t) 6 fj (t) 6 fmaxj (t)
(2)
∀ 1 6 j 6 m, ∀ 0 6 t 6 T
Each contribution fj is constrained by the maximal power output decrease of the associated
machine Mj , which expresses its limited availability. Evolution over time of this maximal
power output, fmaxj (t), is a function of the use of machine Mj , fj (t). Indeed, fmaxj (t)
evolves only if Mj is used, that is, if fj (t) > 0. A first formulation is proposed in Equation (3),
with aj the speed associated to the maximal power output decrease (aj < 0).

fmaxj (t − 1) + aj if fj (t) > 0 ;
fmaxj (t) =
(3)
fmaxj (t − 1)
if fj (t) = 0
∀ 1 6 j 6 m, ∀ 1 6 t 6 T
Equations (2) and (3) being not convex, they can not be used as is within the proposed
convex programming paradigm. A second formulation of the constraints is proposed in set
of equations (4), which details the mathematical program associated to the optimization
problem. This program does not respect the real evolution over time of the maximal power
output that can be delivered by machines, but presents the advantages of being convex
5

and thus consistent with the convex resolution methods proposed in next section. In the
following, the machines behavior follows the simplified model depicted in Figure 2.

!
m
X



min φ σ(·) −
fj (·)
(4a)




j=1




 s.t. fmaxj (t) 6 fmaxj (t − 1)
+µ · aj · (fj (t − 1))υ ∀ 1 6 j 6 m,
(4b)




∀ 1 6 t 6 T, with µ and υ ∈ [0, 1]





with 0 6 fj (t) 6 fmaxj (t)
(4c)



∀ 1 6 j 6 m, ∀ 0 6 t 6 T
Pj (t)
P maxj (0)

P maxj (0)
100%
fmaxj (t) = fmaxj (t − 1) + aj

Decrease
of P maxj
P minj

15%

0

RULmaxj

time

Figure 2: Simplified convex evolution of available power outputs for a machine Mj

3. Resolution based on successive projections onto the sets of constraints
The first proposed resolution approach makes use of convex projections to cope with the
scheduling problem. The general idea is to perform successive projections onto each set of
constraints previously defined in the convex program (4) in order to generate a sequence
of points that converges to the solution of the considered optimization problem. The three
projections are first detailed. The global projection scheme is then described.
3.1. Projections
Each projection allows to comply with one of the three main constraints defined by
Equations (4a), (4b) and (4c).
3.1.1. Reaching of the demand
First constraint corresponds to the reaching of the demand σ(t) at each time t (0 6 t 6
T ). The decision horizon is split into time intervals of length ∆T , with 1 6 ∆T 6 T . For
each time interval, delimited by the times tstart and tend such that tstart 6 t 6 tend, the
projection of the solution F on the demand σ is performed following the strategy detailed
hereafter:
6

Machine selection: machines Mj (1 6 j 6 m) are sorted following an ascending order
of the difference (fmaxj (t) − fj (t)). This sorting favors machines whose contributions
before the projection are the closest from their maximal reachable output fmaxj . This
allows to use first machines already started and to avoid stop-and-start of fuel cells,
which has been shown to lead to premature aging [4].
With distfmax(j) = max(fmaxj (tend) − fj (tend), 0) being the difference between the
current contribution of machine Mj and its maximal reachable contribution, the subscript of the machine Mj 0 whose contribution has to be modified first is determined by
Equation (5).
j0 =
argmin
distfmax(j)
(5)
16j6m|distfmax(j)>0

Projection method: the contribution of the selected machine Mj 0 is increased on each
time interval of length
Pm ∆T just enough to reach the demand σ. The increase being
inc = σ(tend) − j=1 fj (tend), the contribution of machine Mj 0 after the projection
is defined by Equation (6).
fj 0 (t) = min (fj 0 (tend) + inc, fmaxj 0 (tend))
∀ tstart 6 t 6 tend

(6)

The selection of a machine and the projection of its contribution on the demand are
performed on each time interval of length ∆T , until the demand is reached or until there is
no machine available anymore.
3.1.2. Evolution of the maximal output
Each maximal output reachable fmaxj is updated as a function of the contribution of
the corresponding machine, fj , following the relation defined in Equation (7).
fmaxj (t) = fmaxj (t − 1) + µ · aj · (fj (t − 1))υ
∀ 1 6 j 6 m, ∀ 1 6 t 6 T, with µ and υ ∈ [0, 1]

(7)

3.1.3. Respect of the maximal output
A projection of each machine contribution fj onto the corresponding fmaxj can be necessary to fix some possible overrun and to comply with the updated maximal reachable
output. This projection simply limits each contribution fj (t) to the corresponding fmaxj (t)
defined by the application of the previous projection.
fj (t) = fmaxj (t) if fj (t) > fmaxj (t)
∀ 1 6 j 6 m, ∀ 0 6 t 6 T

(8)

3.2. Resolution algorithm based on the projections
The proposed resolution scheme performs successive projections on the three main constraints, following the strategies previously defined. The positivity of each component of the
7

solution vector F is ensured through a positive initialization of each machine contribution
(0)
at each time (fj (t) > 0 ∀ 1 6 j 6 m, ∀ 0 6 t 6 T ).
The evolution of the vectors fmaxj being determined as a function of each fj , a combined
determination of these two elements is required and several successive launches of all the
projections are necessary to optimize the solution. The definition of a schedule that maximizes the production horizon while complying with all the constraints requires then several
iterations of the sequence composed of the three successive projections. This sequence is
iterated until a stopping condition is verified. The number of iterations necessary to obtain
a satisfying solution depends indeed strongly on the configuration of the considered problem, that is, on the number of machines, on the shape of the maximal reachable outputs
fmaxj and on the demand σ. The optimization process of the solution F is then stopped
when it does not change significantly anymore, that is, as soon as the difference between
two successive values of F (|F − F prev |) gets below a certain threshold .
Solutions obtained with this first resolution method are piecewise defined and the reached
production horizons depend strongly on the initialization of the solution vector F and on
the projection strategy. This first resolution approach does then not fully comply with the
resolution paradigm proposed in this paper, which aims to build the solutions globally on the
full production horizon. It allows however to obtain solutions very quickly, which can serve
as a point of comparison with the solutions obtained with the resolution method proposed
in next section.
4. Resolution based on a smooth penalization approach using the Mirror Prox
The second resolution method is based on a smooth penalisation approach. A Mirror
Prox method is proposed to be used to cope with the problem of minimizing the objective
function detailed previously in Equation (4a).
4.1. Principle of the Mirror-Prox
The Mirror Prox algorithm, developed by Nemirovski [11], is a variant of the Mirror
Descent algorithm, which has first been proposed by Nemirovsky and Yudin [12] for convex
programming. It has been extensively studied recently and several relationships have been
discovered between the Mirror Descent scheme and Bregman-proximal methods.
Both approaches are based on the resolution of a primal-dual saddle point problem, which
allows to take constraints into account. The purpose is to minimize a smooth convex function
under constraints. The Mirror Descent algorithm makes use of a gradient descent to find the
minimum of the considered function. The goal is to minimize the local linearization of the
function while not moving too far away from the previous point, with distances measured
via the Bregman divergence of the mirror map. A mirror function allows to transition from
the primal space, where all the constraints of the problem are defined, to the dual space.
The Mirror Prox method applies at each iteration two consecutive steps of Mirror Descent.
A very instructive description of the Mirror Prox algorithm has been proposed by Bubeck
[6]. As illustrated in Figure 3, the Mirror Prox algorithm applies a first time the Mirror
Descent to go from xt to yt+1 and then a second similar step to obtain xt+1 . This second
8

step stats again from xt , but uses the gradient of f evaluated at yt+1 (instead of xt ). These
steps are defined by Equations (9) to (12).
∇θ
∇θ(xt )
−λ∇φ(xt )

−λ∇φ(yt+1 )

∇θ(x0t+1 )
Gradient steps
Eq. (9) and (11)
0
)
Rn ∇θ(yt+1

x0t+1

(∇θ)−1
Dual space

xt
xt+1
yt+1

Projection steps
0
(10) and (12)
yEq.
t+1
D
Primal space

Figure 3: Operating principle of the Mirror Prox method, based on [6]

0
∇θ(yt+1
) = ∇θ(xt ) − λ∇φ(xt )

with λ > 0

0
0
yt+1 = PC (∇θ(yt+1
)) = argmin Dθ (x, yt+1
)

(9)
(10)

x∈C∩D

∇θ(x0t+1 ) = ∇θ(xt ) − λ∇φ(yt+1 )

with λ > 0

xt+1 = PC (∇θ(x0t+1 )) = argmin Dθ (x, x0t+1 )

(11)
(12)

x∈C∩D

4.2. Resolution algorithm based on the Mirror Prox
The proposed resolution algorithm makes use of the Mirror Prox scheme previously described. The resolution of the primal-dual problem being highly sensitive to the variation of
the different parameters associated to the model, the proposed formulation incorporates the
different constraints of the application directly into the objective function. The considered
objective function φ(F, fmax, σ), detailed in Equation (13), aims then at satisfying both the
reaching of the demand and the respect of the evolution of the maximal outputs fmaxj (t),
∀ 1 6 j 6 m and ∀ 0 6 t 6 T .
φ(F, fmax, σ) =λdem hdem (F, σ) + λslope hslope (F, fmax)
with λdem , λslope > 0

(13)

The first component of this objective function, hdem (F, σ), detailed in Equation (14),
aims at satisfying the demand σ(t) at each time t. We did not choose the L2 loss because
asymetric loss are slightly more relevant in our context. In fact, finding a solution which
does not match the demand should be more heavily penalised than a solution which matches
the demand by far. This cost function, based on the objective function previously defined in
Equation (4a), allows to heavily penalize F when the demand becomes larger than the total
output provided by the set of machines and conversely to penalize it only slightly when the
production is larger than the demand. Second part of the objective function, hslope (F, fmax),
is used to control the evolution of each fmaxj as a function of each associated contribution
9

fj (see equation (15)). This corresponds to the penalization associated to the constraint
previously defined by Equation (4b).

X

T
m
X
1
hdem (F, σ) =
exp − γ
fj (t) − σ(t)
T +1
t=0
j=1
with γ > 0

hslope (F, fmax) =

T X
m
X

(14)

 
exp δ fmaxj (t) − fmaxj (t − 1)

t=1 j=1
υ0

0



− µ aj (fj (t − 1))

(15)

with δ > 0, µ0 > 0, υ 0 > 1
The mirror function θ considered in the Mirror Prox formula is defined on Rm(T +1) by
Equation (16) and its gradient, which is used in the Mirror Descent steps, is expressed in
Equation (17).
θ(F ) =

T X
m
X

F ln(F )

(16)

t=0 j=1

∇θ(F ) = ln(F ) + 1

(17)

The Mirror Prox used in the proposed resolution algorithm is then of the following form
(Equations (18) to (21)):
∇θ(F (l+1) ) =∇θ(F (l) ) − λ∇F φ(F (l) , fmax, σ)
with λ > 0

(18)


Fint = exp ∇θ(F (l+1) ) − 1

(19)

∇θ(F (l+1) ) =∇θ(F (l) ) − λ∇Fint φ(F (l) , fmax, σ)
with λ > 0

(20)


F (l+1) = exp ∇θ(F (l+1) ) − 1

(21)

with the gradient of φ defined by Equations (22) to (24).
∇F φ(F (l) , fmax, σ) =λdem ∇F hdem (F, σ)
+ λslope ∇F hslope (F, fmax)
with λdem , λslope > 0
10

(22)

T
X

1
·
T +1
t=0

X

m
exp − γ
fj (t) − σ(t)

∇F hdem (F, σ) = − γ

(23)

j=1

with γ > 0

0 0

∇F hslope (F, fmax) = −δµ υ

T X
m
X

0

aj · fj (t − 1)υ −1 ·

t=0 j=1

 

0
υ0
exp δ fmaxj (t) − fmaxj (t − 1) − µ aj (fj (t − 1))
with δ > 0, µ0 > 0, υ 0 > 1

(24)

In order to accelerate the convergence of the method, the distance of the solution to
its fictional projection onto the demand σ(t), (F − F proj ), is proposed to be added to each
gradient step performed by the Mirror Descent steps. This allows to accelerate the evolution
of the solution towards a state which complies at best with the different constraints. The
new formulation of the Mirror Prox is then detailed in Equations (25) to (28), with the
vector F proj determined by a projection of the solution vector F onto the demand σ(t), with
gF proj = wgrad (F (l) − F proj ) and wgrad > 0.
∇θ(F (l+1) ) =∇θ(F (l) )
− λ ∇F φ(F (l) , fmax, σ) + gF proj
with λ > 0




Fint = exp ∇θ(F (l+1) ) − 1

(25)

(26)

∇θ(F (l+1) ) =∇θ(F (l) )
− λ ∇Fint φ(F (l) , fmax, σ) + gF proj
with λ > 0

F (l+1) = exp ∇θ(F (l+1) ) − 1


(27)
(28)

The use of this additional gradient step allows furthermore to define a stopping condition
for the global process, which allows to avoid the definition of a number of iterations adapted
to the method and to the optimization problem parameters. The resolution process is
stopped when the differences between all the components of the solution vector F and the
corresponding components of the vector F proj are all lower than a certain threshold. In
11

other words, the resolution process is stopped when the reaching of the demand can not
be improved anymore. After each step of the Mirror Prox, several projections onto the
constraints are moreover performed in order to guarantee the respect of these constraints
and to improve the convergence speed of the solution. An update of the maximal power
outputs reachable, fmaxj , as a function of each vector fj , allows to accelerate their evolution
even if their progression is handled by the Mirror Prox. A projection of the machines
contributions fj onto the corresponding fmaxj allows finally to comply with the maximal
power output limitation.
5. Simulation results
5.1. Problem generation
Random problem configurations have been generated using a simulator and configured
with many parameters including the number of machines in the considered multi-stack fuel
cell system, m, and intrinsic fuel cell characteristics. The latter have been defined on the
basis of fuel cell manufacturer specifications and considering a maximal lifetime RULmaxj =
1500 hours ±20% for each machine Mj (1 6 j 6 m). Each RULmaxj value is drawn following
a uniform distribution between 1200 and 1800 hours. Power output values are determined
in the same way, with P maxj (0) = 500 W ± 5% and P minj = 0.15 · P maxj (0) for each
machine Mj .
For the results presented hereafter, the power demand has been assumed to be constant
during the whole scheduling horizon: σ(t) = σ ∀ 0 6 t 6 T . Without any loss of generality,
only one demand value has then been associated to each problem configuration, but many
demands corresponding to different configurations have been tested. Many loads α have been
defined such that σ = α · P nomtot , with P nomtot the nominal total powerPoutput reachable
with the considered fuel cells system and 30% 6 α 6 90%. P nomtot = m
j=1 P nomj , with
P nomj = 0.75 · P maxj (0) the power output recommended by fuel cell manufacturers for a
nominal use of fuel cells. In the following figure, results are represented as a function of the
load α.
5.2. Resolution methods configuration
Initial values of each solution vector F has first been set to zero: fj (t) = 0 ∀ 1 6 j 6 m,
∀ 0 6 t 6 T . Quality of solutions from the point of view of the reached production horizon
globally increases with the number of iterations and stabilizes starting from a certain value.
For each resolution method proposed earlier, several iterations of the associated process are
then performed to optimize this solution. The global process is stopped when the variation
of the solution vector is not significant anymore, that is when (F − F proj ) <  for the
Mirror-Prox-based algorithm (resp. (F − F prev ) <  for the algorithm based on the
successive projections), with the threshold  defined as a function of the demand value as
follows:  = 0.1 · σ.
Tuning of the different parameters involved in the two resolution methods allows to
comply with the constraints and to adapt the shape of each fmaxj . For the evolution of
the maximal power output fmax, values have been defined as follows: µ = 0.2, υ = 0.3.
12

Values for the parameters used in the Mirror-Prox algorithm are the following: λ = 5 · 10−5 ,
λdem = λslope = 100, γ = 100, δ = 100, µ0 = 1, υ 0 = 1.
5.3. Post-processing
The main constraint of the optimization problem is the reaching of the power output demand σ(t). This constraint being tackled through the minimization of an objective function,
solutions may contain time periods during which this demand is not reached. But, with the
two resolution methods proposed in previous section, solutions are built so that the time
periods for which the power demand σ(t) is reached are gathered at the beginning of the
schedules. This is consistent with the objective to maximize the production horizon of the
set of machines. This behavior is linked with the shape of the functions fmaxj (t) representing for each machine Mj the evolution over time of the maximal power output reachable,
which limits the contribution of each machine. These functions being strictly decreasing
with the use of machines, it is in fact more likely to reach the demand at the beginning of
the scheduling process than after some time. As already mentioned, resolution algorithms
detailed previously can then be applied on overestimated horizons T , named decision horizons. The production horizon of each solution, H, is simply the maximal time during which
all the constraints are strictly satisfied. In practice, the production horizon corresponds to
the time during which the demand σ(t) is reached.
5.4. Results
Efficiency of the proposed commitment strategies defined in Sections 3 and 4 is assessed
through the comparison of reached production horizons to a theoretical upper bound. Considering a constant demand σ and a set of fuel cell stacks, this upper bound, denoted UB
and defined in Equation (29), corresponds to the theoretical maximal time during which the
demand can be reached.

P

 m
 0.6 · P maxj (0) · RULj (P minj ) 

 j=1

(29)
UB = 


σ
Figure 4 shows the production horizons obtained with the considered strategies normalized with the upper bound, when considering a set of 25 fuel cell stacks. This upper bound
being never reachable, results are actually better than showed. One can see that the first
strategy, performing successive convex projections, allows to reach a mean relative horizon
of around 39% of the upper bound UB and does not exceed 51, 6%. The resolution method
based on the Mirror Prox gives better results. It allows indeed to reach a mean relative
horizon of around 64.3% of UB and 74, 7% for the best case.
The two resolution approaches differentiate also themselves by the shape of the solutions
and by their computation time. Figures 5 and 6 show schedules obtained respectively with
the algorithm based on the successive projections and with the one based on the Mirror
Prox, when considering m = 3 machines. Evolution of each machine contribution and
13

Normalized horizon (%)

100

75

Resolution method
Projections
50

Mirror Prox
UB

25

0
30

40

50

60

70

80

90

load (%)
Figure 4: Normalized horizon – m = 25 machines

of the associated maximal power output reachable, fmaxj , is also shown in these figures.
One can see in Figure 6 that the algorithm based on the Mirror Prox allows to reach a
better production horizon (H = 1745 periods of time) than the one performing successive
projections onto the sets of constraints (H = 1147). The Mirror Prox-based algorithm
defines also a smooth use of machines, which complies with a continuous use, without sudden
change of output such as those that can be seen in the solutions obtained with the successive
projections (see Figure 5).
The method performing successive projections onto the sets of constraints is not efficient
in terms of reached production horizons, but allows to obtain solutions in very limited time
(less than 9 seconds1 for all the scenarios tested when considering 25 machines). The time
needed to obtain satisfying solutions with the Mirror-Prox-based algorithm is longer (5.6
minutes on average), but does not exceed 17 minutes1 for 25 machines (see Figure 7) and
35 minutes1 for 100 machines, for all the power demands.
6. Conclusion
A management of fuel cell systems has been proposed in a PHM framework. Decision
coming within the scope of Prognostic Decision Making has been addressed considering
longer time frames than those proposed so far in the literature on fuel cells. The use of
1

Simulations have been made using Matlab (Parameters:
CPU@3.30GHz×4, 15.6 Gio, 64 bits)

14

Processeur Intelr CoreT M i5-3550

convex programming has been proposed to cope with the scheduling problem of multi-stack
fuel cell systems under service constraint. A mathematical formulation of the problem has
been proposed as well as two different convex resolution methods performing a minimization
of the objective function under constraints. First one is based on successive projections onto
the sets of constraints and second one on the Mirror-Prox algorithm.
All the fuel cell properties are not observed by the solutions obtained with the proposed
approaches, but this first study is promising. It shows indeed that a global resolution on the
scale of the whole production horizon can be used to define the commitment of machines
over time with the production horizon maximization as objective.
Acknowledgements: This work has been supported by the Labex ACTION project
(contract “ANR-11-LABX-0001-01”)
References
[1] Balaban E and Alonso JJ (2012) An approach to prognostic decision making in the aerospace domain.
In: Annual Conference of the Prognostics and Health Management Society.
[2] Bauschke HH and Borwein JM (1996) On projection algorithms for solving convex feasibility problems.
Society for Industrial and Applied Mathematics (SIAM) Review 38 (3): 367–426.
[3] Beck A and Teboulle M (2003) Mirror descent and nonlinear projected subgradient methods for convex
optimization. Operations Research Letters 31: 267–175.
[4] Borup R, Meyers J, Pivovar B, Inaba M, Ota K, Ogumi Z and et al (2007) Scientific aspects of polymer
electrolyte fuel cell durability and degradation. Chemical reviews 107(10): 3904–3951.
[5] Bosco AD and Fronk MH (2000) Fuel cell flooding detection and correction. US Patent nb. 6.103.409.
[6] Bubeck S (2014) Theory of Convex Optimization for Machine Learning. arXiv preprint arXiv:1405.4980.
[7] Herr N, Nicod JM and Varnier C (2014) Prognostic decision making to extend a platform useful life
under service constraint. In: Proceedings of IEEE International Conference on Prognostics and Health
Management (PHM). Spokane, Washington.
[8] Herr N, Nicod JM and Varnier C (2014) Prognostics-based scheduling in a distributed platform: Model,
complexity and resolution. In: 2014 IEEE International Conference on Automation Science and Engineering (CASE), Taipei, Taiwan. pp. 1054–1059.
[9] Herr N, Nicod JM, Varnier C, Jardin L, Sorrentino A, Hissel D and Péra MC (2017) Decision process
to manage useful life of multi-stacks fuel cell systems under service constraint. Renewable Energy 105C:
590–600. DOI:10.1016/j.renene.2017.01.001.
[10] Jouin M, Gouriveau R, Hissel D, Péra MC and Zerhouni N (2013) Prognostics and health management
of PEMFC - State of the art and remaining challenges. Int. Journal of Hydrogen Energy 38: 15307–
15317.
[11] Nemirovski AS (2004) Prox-method with rate of convergence O(1/t) for variational inequalities with
lipschitz continuous monotone operators and smooth convex-concave saddle point problems. SIAM
Journal on Optimization 15 (1): 229–251.
[12] Nemirovsky AS and Yudin DB (1983) Problem complexity and method efficiency in optimization. Wiley
- Interscience Series in Discrete Mathematics.
[13] Saha B, Quach CC and Goebel K (2012) Optimizing battery life for electric UAVs using a bayesian
framework. In: Proc. of IEEE Aerospace Conference, Big Sky, MT, USA.
[14] Wells B and Parr RK (2004) Fuel cell system method, apparatus and scheduling. US Patent nb. 0185316
A1.

15

600

600

fmax1(t)

fmax2(t)

400

200

0
0

M2 contribution (f2(t))
Output P(t)

Output P(t)

M1 contribution (f1(t))

1000
2000
Time (hours)
600

400

200

0
0

3000

1000
2000
Time (hours)

3000

fmax (t)
3

Output P(t)

M3 contribution (f3(t))
400

200

0
0

1000
2000
Time (hours)

3000

(a) Machines contributions
600

Output P(t)

500

H = 1147

400

Demand σ
M3 contribution (f3(t))

300

M2 contribution (f3(t)+f2(t))
M1 contribution (f1(t)+f2(t)+f3(t))

200
100
0
0

500

1000

1500
Time (hours)

2000

2500

3000

(b) Schedule
Figure 5: Solution obtained with the successive projections – m = 3 machines, σ = 0.4 · P nomtot

16

600

600

fmax1(t)

fmax2(t)

400

200

0
0

M2 contribution (f2(t))
Output P(t)

Output P(t)

M1 contribution (f1(t))

1000
2000
Time (hours)
600

400

200

0
0

3000

1000
2000
Time (hours)

3000

fmax (t)
3

Output P(t)

M3 contribution (f3(t))
400

200

0
0

1000
2000
Time (hours)

3000

(a) Machines contributions
500

H = 1745

Output P(t)

400
Demand σ
M3 contribution (f3(t))

300

M2 contribution (f3(t)+f2(t))
M1 contribution (f1(t)+f2(t)+f3(t))

200

100

0
0

500

1000

1500
Time (hours)

2000

2500

3000

(b) Schedule
Figure 6: Solution obtained with the Mirror-Prox-based algorithm – m = 3 machines, σ = 0.4 · P nomtot

17

Computation time (s.)

750

Resolution method

500

Projections
Mirror Prox

250

0
30

40

50

60

70

80

90

load (%)
Figure 7: Computation times – m = 25 machines

18

