Determination of the Equivalent Consumption in Hybrid Electric Vehicles in the State-Constrained Case

— In the optimization of power management of hybrid electric vehicles, the equivalent consumption factor is often used. This parameter represents a way of penalizing the use of power from the batteries, taking into account the fuel consumption that such use eventually hides. If the problem of determining the power split between the energy sources of the vehicle that minimizes fuel consumption is stated as a non linear constrained optimal control problem, and is solved using Pontryagin Maximum Principle (PMP), the equivalent consumption factor may be computed from the adjoint state. Following this approach we compute the trajectory of the adjoint state in the case where state constraints are taken into account. The optimality conditions from PMP are a Boundary Value Problem (BVP), which is solved numerically using a code named PASVA4. Numerical examples are compared with dynamic programming solutions of the same problem. It is found that the adjoint state is continuous and its trajectory is described. The approach may be generalized to similar optimal control problems.


INTRODUCTION
In the operation of hybrid electric vehicles, in order to actually achieve fuel savings, a control system that determines the power split among the energy sources according to the driver requirements must be defined.This control, usually electronically implemented, is called "supervisory control''.Each one of the vehicle devices has in turn its own controller that obeys the supervisory commands at a lower level.
The objective of the supervisory control is to minimize fuel consumption along a given mission.In the case of public transport, where the routes are repetitive or in vehicles equipped with positioning systems, it is possible to know in advance the future power requirements.In the general case, there is the drawback that the future power requirements are not known.
Thus, many real time supervisory control strategies are based in the so called "Equivalent consumption" [1][2][3].The idea is to find a weighing factor that penalizes the instantaneous energy consumption from the batteries taking into account that, as the batteries will be recharged by means of the engine, any use of energy from the batteries eventually hides a fuel consumption.If this weighing factor is known at each time, the problem of minimizing the energy consumption in a duty cycle, a global problem, can be transformed in a local minimization problem and thus, can be solved in real time.Clearly this weighing factor varies with the required velocity cycle and the relative efficiency of the engine and electrical paths of the vehicle.Many researchers have studied the form of determining this parameter [1][2][3][4][5].
In previous work [6], we have formulated the supervisory control problem of hybrid vehicles as a non-linear optimal control problem with control and state constraints.It have been also stated the optimality conditions given by the PMP [7][8][9][10][11][12][13]. Within this approach, it can be seen that the equivalent consumption factor is related to the adjoint or co-state of the optimal control problem [1].So that, by solving the optimality conditions and obtaining the trajectory of the adjoint state the instantaneous value of the equivalent consumption factor could be found.This is then our aim.
It must be pointed out that the optimality conditions are differential-algebraic equations and thus their numerical solution may not be straightforward [14,15].Sometimes, the algebraic equations may be solved as a function of the dynamic variables and the problem may be reduced to solve an ordinary differential equations system with boundary values.If there are constraints, discontinuities or switchings in the right hand side of the differential equations may appear.Moreover, they may appear in times that are not known in advance but depend on the very solution [6].
PASVA4 is a boundary-value non-linear ordinary differential equations system solver which is capable of handling discontinuities in the equations and in the solutions and solving simultaneously for algebraic equations and unknown parameters [16,17].These capabilities make it useful to solve this optimal control problem.
Pontryagin Maximum Principle (PMP) has been used before to address similar supervisory control problems in [5,11,[18][19][20], etc.In these articles, in order to consider the state constraint, the optimality conditions are formulated following the "indirect adjoining approach''.The adjoint state defined in the frame of this approach may experience jump discontinuities and the optimality conditions are not sufficient to determine these jumps.Then, these authors determine the values for the adjoint state by a dichotomic search.In [21], we presented a numerical form of computing the adjoint state in the case where the state constraint is not reached.In this work, we extend those results to the case where the state constraint becomes actually active.We found that if other forms of adjoining the state constraint are used, the optimality conditions allow the computation of the complete trajectory of the adjoint state as a continuous variable.
In [22], an algorithm is proposed which successively solves problems where the state constraint is removed.Authors show, for two simplified problems, the convergence of these successive solutions to the original one, and then extend results for a model of a hybrid vehicle.This approach has the advantage that it does not need to know in advance the number of times where the state constraints are active.Nevertheless, the authors recognize that, "the range of applicability of the method remains to be clarified''.We will return to these points in Section 3.
Another approach to a similar problem has been recently presented in [23], where unconstrained problems are iteratively solved, by means of a shooting algorithm.The authors perform consecutive subdivisions in the intervals where the state constraints are exceeded, forcing the boundary conditions in the subintervals to fulfil the constraints in such a way that convergence to the state constrained trajectory is obtained.This approach neither need to know a priori the number of binding intervals and is more general than ours in regard to the assumptions made on the functions involved (Sect.3).We believe that it would be meaningful to analyse if the application to our particular problem lead to similar results.
We point out that these two last approaches use shooting methods to solve the involved BVP.Instead we use a global finite difference method, known to be more stable for general problems.
In [24], it is pointed out that the approach based on the PMP allows an unified treatment of supervisory control problems of increasing difficulty.Different variants of the basic supervisory control problem are discussed, namely the addition of new terms to the objective functional (as pollutant emission rates, the integral of current flow in and out the batteries, which is related with battery aging), modifications to boundary conditions (as condition of battery depletion for plug-in hybrid electric vehicles), addition of state variables (as engine temperature, engine and catalyst temperature, emmissions in Diesel engines, battery temperature) and addition of new state constraints.This last point is particularly an important issue in the case of hybrid electric vehicles equipped with batteries and ultracapacitors since energy in the ultracapacitors is a new state variable that, because of the high power that ultracapacitors may deliver, will frequently reach its lower limit.Hence, the approach proposed in this article will contribute to the solution of all these problems.

MODEL AND PROBLEM STATEMENT
The model and optimal control problem has been described in detail in [6].Here, we shall only recall that we used a simplified model of the vehicle that only takes into account the energy sources and the power flows among them [25,26].The control action uðtÞ is the power delivered by the engine.The sum of the power delivered by both sources must always equal the power required by the driver which is called rðtÞ.Hence, the power supplied by the batteries is rðtÞ À uðtÞ.Concerning fuel consumption, it must be taken into account that, to compensate for the losses that occur at each one of the energy conversion processes involved in the path from the fuel tank to the wheels (chemical to mechanical, mechanical to electrical, etc), for the source to effectively provide a power uðtÞ, it must actually deliver a greater amount.This effect will be represented by a function f C that depends on u itself.Fuel consumption is then represented by the integral of the total delivered power, i.e., R T 0 f C ðuÞdt.A similar situation occurs in the electrical path, and therefore the energy within the batteries at each time is given by: where x 0 is the initial energy, f B is the function that represents the compensation for losses in the electrical path and xðtÞ is the energy in the batteries at time t, which is considered the state variable.We impose that the final energy xðT Þ be equal to the initial, which means a charge sustaining operation of the vehicle.The functions f C and f B are, in practice, determined by tests performed on the vehicle at different operation points.In this work, we will assume that these functions are determined by approximation of the experimental points in such a way that they are continuously differentiable.We also assume that f B ð0Þ ¼ 0 and that f B ðsÞ !s for all possible value of power from the batteries path.
There are constraints in the power flows because of the physical limitations of the devices involved, and a lower bound for the state variable, which says that the batteries cannot be discharged below certain minimal value x min , in order to avoid any damage.Then, the problem is the following: Find a piecewise continuous function u such that maximizes: subject to: x min xðtÞ ð 7Þ where K max , K min , u max , x min are vehicle parameters and x 0 is the initial energy in the batteries.

DIRECT ADJOINING APPROACH OPTIMALITY CONDITIONS
In order to solve this problem using maximum principle optimality conditions, we define the Hamiltonian: where k is the adjoint state.In order to take into account the constraints, we define the Lagrangian.Let us define the functions U ðtÞ ¼ maxf0; rðtÞ À K max g and U ¼ minfu max ; rðtÞ À K min g and hðxÞ ¼ x À x min to represent the constraints on control as U ðtÞ u U ðtÞ and the constraint on the state as hðx; tÞ !0.
There are different forms of adjoining the state constraint to the Hamiltonian.The most commonly used is the "indirect adjoining'' form, which consists of adding a multiplier times the first total time derivative of h where the control appears explicitly.In this case, this time derivative is ðr À uÞ.However, in previous work [6] we showed the difficulty to determine, within this approach, the jumps that may occur in the adjoint state when the state enters or leaves the binding interval.This fact is also mentioned by other researchers, who overcome this drawback by determining the jumps by a dichotomic search [5,11,13,18].
To overcome this drawback, here we use the "direct adjoining'' approach as stated in [9].In this form, the optimality conditions lead to a solution where the adjoint state is continuous and that coincides with the solution found using dynamic programming.In this approach, the Lagrangian is: where hðtÞ, hðtÞ and lðtÞ are Lagrange multipliers.
In [9], it is presented a third approach for adjoining the state constraint.It consists of adjoining indirectly but assumming that the adjoint state is continuous.The Lagrangian in this case is: Note that the adjoint state and Lagrange multipliers of the different forms of adjoining are different, so they are denoted with different characters.Nevertheless there is a relation between them.In [9], it is shown that: If the optimality conditions in this approach are stated and solved, the obtained k and l satisfy Equation (11).In addition, recalling that our purpose is the determination of an instantaneous factor that "weighs'' the power from the electric path in order to set it in terms of fuel consumption, it turns up that k in H and k þ l in H represent that factor.Hence, in the next we will state the optimality conditions under the direct adjoining approach.
We assumed that the state constraint becomes active only in a subinterval ½t 1 ; t 2 , interior to the interval of interest.The optimality conditions are then the following [9]: At s ¼ t 1 and s ¼ t 2 : In Equation ( 13), f 0 C and f 0 B indicate respectively the derivatives of f C and f B with respect to their unique variable, i.e., f 0 C ðsÞ ¼ df C ðsÞ ds and f 0 B ðsÞ ¼ df B ðsÞ ds .

SOLUTION OF THE OPTIMALITY CONDITIONS
Beginning at the final time T , assuming xðT Þ 6 ¼ x min and using Equation (20), we obtain c ¼ 0 and then kðT À Þ ¼ b.
In the interval ðt 2 ; T , x > x min and so, using Equation ( 18), l 0. Using Equation ( 14), it results that k is constant and k b.In this interval, uðtÞ and u 6 ¼ U , because of Equation ( 12), u must be solution of: Then u may be expressed in terms of k: If t 2 were known, the value for k, (equivalently b), could be found by solving the two point BVP given by: However, t 2 is not known a priori and an additional equation is needed to obtain it.
In the boundary interval ðt 1 ; t 2 Þ, x x min , and then _ x ¼ 0 and Àf B ðr À uÞ ¼ 0. Because of the form we assumed for f B , the latter implies that u r in that interval.In other words, the energy in the batteries has reached its lower bound and then all the power must be provided by the engine.By Equation (21), for s ¼ t 2 we have that: Equation ( 26) relates kðt þ 2 Þ to t 2 and then is the additional equation that allows to determine t 2 by solving simultaneously with Equation (24).
From Equation ( 13), once known k and u, h and h may be computed.In ðt 1 ; t 2 Þ, as the control action is obliged to be u r, the value for k losses its significance as an equivalent consumption factor.Nevertheless, let us analyze the optimal adjoint state trajectory in this interval.In the case that U < r ¼ u < U , h ¼ h ¼ 0 and using Equation ( 13): Then, Using Equation ( 25), gðt 2 Þ may be obtained.Note that: is always a solution for Equation (26).Replacing this value in Equation ( 22) for Wherever u ¼ r ¼ U , it must be r ¼ 0 (r cannot be equal to r À K max ).Then using Equation ( 13): Analogously, if u ¼ r ¼ U , it must be r ¼ u max (r cannot be equal to r À K min ).Then using Equation ( 13): and It remains to determine the values for h and h in these limit cases which are the only situations where the adjoint state could have a jump.
From Equation ( 14), in ðt 1 ; t 2 Þ: where f 00 C ðsÞ ¼ d 2 f C ðsÞ ds 2 .At t 1 , an analogous situation to that in t 2 occurs.From the conditions (21), it results: In the interval ½0; t 1 Þ, once again x > x min and then l 0 which means that k kðt À 1 Þ.The control u satisfies uðtÞ ¼ arg max u2X H ¼ uðt; kðt À 1 ÞÞ.Replacing in Equation (36), yields: Gathering this equation and the differential system: we are able to find the solution by means of PASVA4, obtaining simultaneously the unknown value for t 1 .
Note again: is always a solution for Equation (37).If and from Equation (13): , k is continuous at t 1 .Table 1 summarizes the above results.Note that the adjoint state k takes, in the interval ½0; t 1 a constant value, varies in the binding interval as a function of the required power, and it is again constant in ½t 2 ; T .Because of the continuity, the constant values outside the binding interval are those of the entry and exit times.Note also that we have used the fact that f 0 C and f 0 B exist, which implies that if the efficiency functions are given by a table of experimental values they must be approximated with at least that smoothness.We have also used that at s ¼ 0, f 0 B ðsÞ is different from 0. Another key point is the assumption that the function f B is zero at the origin, which is sustainable since there will be no additional losses if there is no power supply from the batteries.
As it has been said, the correctness of this solution may be checked by solving the optimality conditions from the "indirect adjoining with continuous adjoint state'' approach.That derivation is presented in the Appendix.
If a similar procedure is derived for the model presented in [22], similarities are found between the expressions for the trajectories previous to the binding interval.Unfortunately, the behaviour of the adjoint state in the whole interval is not explicitly shown.In [23], concerning the trajectory of the adjoint state, authors start from the assumption that it presents jumps in the boundary interval and then, they compute these jumps according to a more general theorem than the one we used, that is also presented in [9].It would be very enlightening to apply their algorithm to our particular model in order to see if the final discrete sequence of the of the adjoint state trajectory with the corresponding jumps approximates the continuous one we found, by using the optimality conditions of the direct adjoining and the indirect adjoining with continuous adjoint theorems [9].We will address this in the future, even though it is difficult to match the many choices that must be decided during the implementation and that do not appear explicitly in their text.

NUMERICAL SET UP
The values of the adjoint state k outside the binding interval are obtained numerically as the solutions of two free boundary problems through the use of PASVA4.This tool is a nonlinear two-point boundary-value solver for first order systems of ordinary differential equations.It is capable of efficiently handling internal jump discontinuities, even when they are at locations that depend on the solution itself, and solving simultaneously for parameters and additional algebraic equations [16,17].It is based on finite differences, i.e., the trapezoidal rule, on a general mesh.The order is enhanced by deferred corrections.It has also a procedure for automatically choosing an appropiate mesh, which is based on the equidistribution of an estimate of the local truncation error.Hence, it is a variable order, variable step method.To handle jump discontinuities, each interval of continuity is mapped into a unit interval where a transformed system is solved.Special precautions are taken for the deferred corrections procedure not to involve mesh points at both sides of the discontinuities and to locate the discontinuity conditions such as not to disturb and to preserve the structure of the linear equations solver at each step of the Newton's method used to solve the linearized discrete equations.The capability of solving simultaneously for unknown parameters is used in this case by considering the free boundary as an unknown parameter.This software tool has been used succesfully for a long time for different problems.
If the control function uðt; kÞ cannot be expressed as an explicit function of k, a numerical minimization of the Hamiltonian using the current values for kðtÞ and xðtÞ may be performed, at each time t in order to obtain uðtÞ.Then, the free boundary problems in Equations ( 24) and (38) may be solved.Note that this is, actually, an algebraicdifferential system and then it may present severe difficulties to be solved [14].Nevertheless, in this case, the procedure of consecutively solving first the optimization problem and then the BVP led to the solution.
To test this procedure, we performed several tests using MINOS 5.4, which is a well-known FORTRAN-based software tool for non-linear programming problems [27,28].
_ r Actually, it has been designed for large scale non-linear constrained optimization problems.In this case, most of the abilities offered by MINOS are not used, but only its capability to solve minimization problems in bounded intervals (i.e., within the set X ¼ fu=U ðtÞ u U ðtÞg).Our choice was only based in our familiarity with the tool, achieved in previous work [29].MINOS may be used as FORTRAN subroutine (MINOSS).It is called by the routine of PASVA4 that performs the evaluation of the right hand side of the differential equation.The only requirement is that the Hamiltonian be convex.If f C and f B are given by tables, we need to assume that they may be interpolated or approximated in such a way that this condition is fulfilled.

NUMERICAL EXAMPLES
In this section, we will present an example of the analytical results presented in the previous section.In addition we will present the solution of the same problem obtained using dynamic programming.The dynamic programming algorithm used is described in [30].With this purpose, we will begin by considering the following assumptions: there is only one interval ½t 1 ; t 2 & ½0; T where the state constraint is active; -rðtÞ is a sine wave.
The choice of f C and f B ensures the smothness and convexity of the Hamiltonian.Note that f B ð0Þ ¼ 0 and df B ðsÞ=dðsÞ ¼ 1 at the origin.The required power is continuous, differentiable and represents an interval of acceleration (rðtÞ > 0) followed by an interval of deceleration (rðtÞ < 0).The knowledge of an analytical expression for rðtÞ allows its evaluation and the evaluation of its time derivative at any intermediate t, so avoiding the errors that could be introduced by an interpolation.Table 2 summarizes the data used in the simulations.Figure 1 shows the optimal control function found, i.e. the power that must be effectively supplied by the internal combustion engine.Figure 2 shows the corresponding optimal trajectory of the energy in the batteries.Results obtained using dynamic programming are also shown for comparison.Table 3 contains the values of the adjoint state outside the binding interval, the entry and exit times of the binding interval, the corresponding consumption computed by our algorithm and that computed using dynamic programming.
In this case, it can be seen that from Equation ( 26) and Equation (37 respectively before and after the binding interval, while kðtÞ ¼ 2a c rðtÞ þ b c in the binding interval.The control function u, except in the binding interval or when it reaches the bounds, is an affine function of the required power rðtÞ.Table 4 summarizes the analytic results in this particular case. The BVP solver needs initial guesses as starting points for the algorithm.We proceed by solving first simpler problems.First, we solve the state unconstrained problem under the assumption that there are not losses in the electrical path.So the function f C is the identity function.We found that for this problem the solution can be found by starting from the zero state as starting point [6].This gives us an initial  Optimal control obtained using PASVA4 (full line) and using dynamic programming (circles); the required power (dashed line) and bounds on control UðtÞ and U ðtÞ (dotted lines) are also indicated.
guess for the state and for the adjoint state, that in that case is constant in the whole interval.In a second step, we solve the problem including the function f C as it is, but disregarding the constraints on the state.Finally, from these results, we look for initial guesses for the entry and exit times by determining the interval where the state constraint is violated in the unconstrained solution.
To extend the computations to the case where the required power rðtÞ does not have an analytical expression but is obtained from the inverse longitudinal dynamics of the vehicle applied to a standardized velocity cycle, it is necessary to include two additional algorithms.First, it is necessary to compute numerically the time derivative r 0 ðtÞ, which is used by PASVA4.This was performed using centered differences.Secondly, as PASVA4 fits the time step to control the global error, an interpolation algorithm must be used to evaluate rðtÞ at the new intermediate points used by PASVA4.In this case, we simply used a linear interpolation routine.
Concerning the assumption of the existence of only one interval where the state constraint is active, our proposal to remove that assumption is setting PASVA4 as if there were a binding interval each N fict points in the mesh.As the software finds automatically the extremes of the binding intervals, if for any interval, say ½t i ; t iþ1 , it results t i ¼ t iþ1 , we conclude that this binding interval is fictitious and that, actually, the state constraint did not reach the lower bound.Next, we show preliminary results obtained in such a situation.In this case, the data are the same used for the previous figures, but the bound on the state was decreased (x min ¼ 0) in order to make the state not to reach the limit.The binding interval collapsed at t 1 ¼ t 2 ¼ 11:40 s. Results are presented in TABLE 4 Behavior of state, control and multipliers in each interval, in the case of quadratic f B and f C Optimal state trajectory obtained using PASVA4 (full line) and using dynamic programming (circles) corresponding to the situation in Figure 1.Figures 3 and 4. Clearly, this proposal must be investigated further, since another strategies can be designed but the knowledge of at least a maximum number of binding intervals is needed to be known.This is a limitation that is not present in the algorithms presented in [22] and [23].
Next, we show the results obtained on a standardized velocity cycle, namely, the cycle US-EPA-Highway Fuel Economy Test, scaled properly in order to show a state trajectory with an interior binding interval, as it was assumed in the theoretical derivation of the solution.Figure 5 shows in dashed line the power required by the vehicle with the parameters in Table 5 to follow this velocity cycle, and in full line the optimal control function found and the corresponding one obtained using dynamic programming (circles line).The bounding functions U ðtÞ and U ðtÞ are also indicated (dotted lines).Figure 6 shows the corrresponding state trajectories and Table 6 the values for k, the binding interval entry and exit times and consumption.Power (kW) Optimal control using PASVA (full line) when it is assumed that a binding interval exists, but it is shown to be fictitious.
Required power (dashed line) and bounds on control (dotted lines).Power (kW) Optimal control using PASVA (full line) and using dynamic programming (circles), for the cycle US-EPA-Highway Fuel Economy Test.Required power (dashed line) and bounds on control (dotted lines) are also included.

CONCLUSIONS
We have used PMP to obtain the trajectory of the adjoint state with the purpose of designing supervisory control equivalent consumption minimization strategies.In this work, we studied the case where the state constraint becomes active, which means that the energy in the batteries reaches the lower limit.For the model used, which assumes some smoothness of the "efficiency functions'', we found that the adjoint state is constant outside the interval in which the state constraint becomes active, and its value depends on the entry and exit times to that interval.In the boundary interval, if the power contribution of each source does not reach the limits, the adjoint state varies as a function of the required power.Analysis of different forms of adjoining the state constraint to the Hamiltonian as summarized in [9], was what led to this detailed description.We think that, although several issues still remain to be completed (as the cases of repeated binding intervals, non-smoothness or state dependence of the efficiency function of the electrical path), our results are a novel contribution that will be helpful in the design of supervisory control strategies.Particularly, the solution of this state-constrained problem applies to vehicles equipped with ultracapacitors, where the energy is another state variable, which will probably reach its lower bound repeatedly during a duty cycle, since ultracapacitors are intended to perform quick and deep discharges during peaks of power demand.
Solving the optimality conditions leads to two BVP, each one with a free boundary.The value at this boundary is obtained simultaneously with the solution.This is possible by means of using a software tool as PASVA4, capable of handling problems with these features.
Although our main purpose was to find the behaviour of the adjoint state, as solution of the optimal control problem gives also the optimal control action, the approach may also be used in the design of model predictive control schemes.The execution times are short enough to allow the repeated solution of the optimal control problem if adequate (neighboring) initial conditions are used.
Moreover, as it has been largely described in [24], the Maximum Principle approach is a unified frame where different supervisory control problems may be treated.Hence this approach and the versatility of PASVA4 may be used together to solve problems with: additional optimization criteria as pollutant emissions and battery aging; new state variables as engine and catalyst temperature, pollutant emissions and battery temperature; different boundary conditions as battery depletion for plug-in electric vehicles; additional state and control constraints.
The iterative linking of PASVA4 with non-linear programming software tools may also be used for other control problems where optimality conditions are differentialalgebraic equations such that the control action cannot be obtained from the algebraic equations as an explicit function of the state variables.
Optimal state trajectory obtained using PASVA4 (full line) and obtained using dynamic programming (circles).Here, we present the derivation of the solution when the state constraint is adjoined indirectly but including the assumption that the adjoint state k is continuous [9].We recall the Hamiltonian in this approach: A.

Optimality Conditions
The optimal solution satisfies [9]: At the final time T A.

Solution
We assume again that the state constraint is active only in an interval ðt 1 ; t 2 Þ & ½0; T .We begin by the interval ðt 2 ; T .Using Equation (A10) and (A4): From Equations (A11) and (A9), l 0. From Equations (A2), (A3) and (A5), with the same arguments used before, the function u that minimizes H is: where uðt; kÞ is the solution of Àf 0 C ðuÞ þ kf 0 B ðr À uÞ ¼ 0. From this last expression, using Equation (A4) and the state equation we obtain the system: with boundary conditions: If t 2 were known, the solution of this system would allow obtaining the value for k (or equivalently b).
In ðt 1 ; t 2 Þ, x x min ; then _ x ¼ Àf B ðr À uÞ ¼ 0, and therefore u ¼ r.Again this means that, as the energy in the batteries x is in its lowest value x min , they can not provide more power for traction so it is the combustion engine which must provide full power, and so u must be equal to the required power r.If r > U or r < U , the problem has no solution, since u cannot exceed its lower and upper limits.
As k is continuous, because of Equation (A4), it results that k kðt Using Equation (A3) we obtain: , u is continuous at t 2 and by Equation (A6) l is continuous at t x > 0, since power must go into the electrical storage system for the energy in it to move away from its lowest value.Therefore, as h 1 ðt À 2 Þ ¼ 0, h 1 is discontinuous at t 2 .From Equation (A8), l must also be continuous at t 2 .Hence, Expression (A14) is the aditional algebraic condition that allows the solution of the free boundary system (A12, A13).In ½0; t 1 , l is again constant because of Equation (A9).In addition it is continuous at t 1 , by the same arguments used at t 2 .Therefore: Note that l depends on t 1 which, so far, is not known.As, in this interval, k is also continuous and constant, which, from Equation (A14) is already known.Using Equation (A3), when the maximum is interior to the interval ðU ; U Þ: and, so, u may be obtained from this last expression as a function of l (since k is known) .Hence, the free boundary problem: with the additional algebraic condition (A16) may be solved for x, l and the unknown parameter t 1 .Results are summarized in Table A1.Table A2 exemplifies the entries in Table A1 for the case where f C and f B are quadratics.Solution of _ x ¼ Àf B ðr À uðt; lÞÞ _ l ¼ 0 x(0) = x 0 , x(t 1 ) = x min l(t 1 ) = (f 0 C (r(t 1 )) À f 0 C (r(t 2 )))/f 0 B (0) x min Figure 1 Figure 2 in batt.kWs)

Figure 4 Optimal state trajectory corresponding to Figure 3
Figure 4Optimal state trajectory corresponding to Figure3.