# Costate Estimation of PMP-Based Control Strategy for PHEV Using Legendre Pseudospectral Method

Authors: Hanbing Wei; Yao Chen; Zhiyuan Peng
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Costate value plays a significant role in the application of PMP-based control strategy for PHEV. It is critical for terminal SOC of battery at destination and corresponding equivalent fuel consumption. However, it is not convenient to choose the approximate costate in real driving condition. In the paper, the optimal control problem of PHEV based on PMP has been converted to nonlinear programming problem. By means of KKT condition costate can be approximated as KKT multipliers of NLP divided by the LGL weights. A kind of general costate estimation approach is proposed for predefined driving condition in this way. Dynamic model has been established in Matlab/Simulink in order to prove the effectiveness of the method. Simulation results demonstrate that the method presented in the paper can deduce the closer value of global optimal value than constant initial costate value. This approach can be used for initial costate and jump condition estimation of PMP-based control strategy for PHEV.


---

## Body

## 1. Introduction

A lot of approaches based on PMP (Pontryagin’s minimum principle) theory have been applied to the energy management optimal control for HEVs (hybrid electric vehicle) [1, 2]. Intrinsically, it realizes the minimization of the Hamilton function. For its simple enough instantaneous optimization compared with the Dynamic Programming (DP) method, the PMP-based optimal control strategy has more potential to be implemented in real-time application. Moreover, not only fuel consumption but also emission can been optimized comprehensively by PMP method. In [3], the author integrated the simplified three-way catalytic converter model into the simulation. So the exhaust emission from the catalyst outlet can be included in Hamilton function. In [4], the author extended the PMP-based supervisory control for HEV with a new state reflecting the thermal state of the engine. The power losses due to low engine coolant temperature are therefore taken into account in the optimal control of HEV.The prominent difficulty of successful application of PMP is the appropriate assignment ofinitial costate. It is the derivatives of Lagrangian multiplier for variation of calculus substantially [5]. As practical meanings of the tradeoff between the fuel consumption and battery depletion, the costate value plays a critical role in the control strategy of HEV. How to choose initial costate λ decides the operation mode of PHEV significantly [6]. If λ is chosen as SOCb, PHEV will operate in blended mode, which means the battery will be depleted at the end of the driving cycle. Alternatively, if λ is designated as SOCa, the EMS (energy management system) will try to deplete the battery as soon as possible. Both the engine and the motor will operate cooperatively. Accordingly, the inferior optimal control strategy will be achieved because of the poor efficiency of battery when SOC reach SOCmin. Unfortunately, it is not convenient to designate the exact costate in real driving condition until now despite its significance. Iterative simulation using shooting method for specific driving cycle has been restricted in the offline optimization [7]. In addition, driving patterns prediction in prior [8, 9], traffic condition from the external traffic information like ITS or GPS has been dedicated to development of the control strategy of HEV [10–15].Different from HEV, plug-in hybrid electric vehicle (PHEV) is always willing to discharge the battery at the end of the trip, not only for the Charge Depletion/Charge Sustaining but also for blended control strategy. The battery in PHEV has more opportunity to trigger the bottom of SOC consequently [16]. Therefore, handling inequality constraints in PMP requires more attention because the constraints possibly produce ambiguity, which is frequently called a jump condition. When the jump condition of costate was considered, the optimal control of PHEV turns to be more complicated. In [6], the authors propose the mathematical derivation of an additional condition necessary for the inequality state constraints and deduce the necessity of the different costate for the jump condition before and after. Nevertheless, the paper only contributes to understanding the physical definition of the costate for jump condition, still failing to produce the costate directly.In the paper, the optimal control problem (OCP) of PHEV based on PMP has been converted to nonlinear programming problem (NLP). Costate can be approximated as Lagrange multipliers of NLP divided by the LGL weights. A kind of general costate estimation approach is proposed for specific driving condition in this way. The outline of the paper is organized as follows. Section2 introduces the schematic and characteristic of the vehicle researched. Section 3 describes the control strategy of PMP for PHEV. Section 4 presents the general approach of estimating the costate using Legendre Pseudospectral Method. The method mentioned above has been applied to predefined circumstance and some discussion has been presented in Section 5. Final conclusions are drawn in Section 6.


## 2. Vehicle Model

The configuration investigated in this paper is a two-clutch single-shaft parallel PHEV as shown in Figure1. The release or engagement of the wet-type multiplate clutch decides the operation mode of the vehicle. The one-way clutch is used to assure that the revolution speed of motor never falls behind that of engine when engine cranked. In this way, excessive friction of the clutch can be refrained from the speed discrepancy between motor and engine. Besides that, energy conversion from electric energy to mechanical energy or vice versa is implemented by the control strategy, by which the energy conversion loss decreases the overall efficiency. Therefore, an appropriate optimal control is needed to distribute output between different power sources and reduce the fuel consumption, drivability, and exhaust emission [17]. Further, optimized strategy to deplete state of charge of the battery at destination is an essential issue for PHEV control.Figure 1
Schematic of the PHEV.The relationship between the speed and power delivered by the power source can be expressed as(1)Picet+Pmott=Pwhltωicet=ωmott=ωwhltik,where Pice and  Pmot are the power delivered by the engine and motor. Pwhl are the driving power requested from wheel. ωice,  ωmot, and  ωwhl are the revolution speed of the engine, motor, and wheel, respectively.


### 2.1. Engine Model

The aim of the optimal control strategy of PHEV is to minimize the fuel consumption and deplete the battery capacity over a driving cycle. Fuel consumption of diesel engine is generally modeled as a map for every possible combination of speed and torque. Obviously, the map can only be represented as nonlinear high order function. For simplicity, appropriate Willans line model is usually used to express the function of the engine power and speed. From the characteristic of Willans line model, we can say that the efficiency of the energy conversion device can be modeled by representing the input power as an affine function of the output power and losses.At any given speed, the engine powerPice can be represented as an affine function of the fuel power Pful. The gradient and intercept of each Willans line can be expressed as polynomial functions normally depending on engine speed, by(2)Pfult=e0ωicet+e1ωicetPicet,e0ωicet=e00+e01ωicet+e02ωice2t,e1ωicet=e10+e11ωicet+e12ωice2t,m˙ft=1QLHVe0ωicet+e1ωicetPicet,where e0e1 are the variable coefficient of the fitted Willans model, as shown in Figures 2 and 3. e00,  e01,  e02,  e10,  e11, and e12 are the constants of fitting. The final fitted Willans model is compared with original fuel map in Figure 4.Figure 2
Curve of coefficiente0 with engine speed.Figure 3
Curve of coefficiente1 with engine speed.Figure 4
Comparison of fitted Willans model and original fuel map.


### 2.2. Battery Model

Traditionally, the battery SOC is appropriate for defining the rate change of power energy stored in the battery. The simplified model of battery can be expressed in (3) according to internal resistance model:(3)Ibat=Voc-Voc2-4RPbat2R,where Ibat is the discharge current of battery; Pbat is the power of battery; Voc is the open voltage of battery; R is the internal resistance of battery.The total electric power loss over the battery resistance is(4)PlPbat=RIbat2=Voc-Voc2-4RPbat24R. Assuming the power loss of battery as constant, the power loss can be approximated by formula of expanded Taylor series around Pbat=0, yielding(5)PlPbat≈Pl0+∂Pl∂Pbat0+∂2Pl∂Pbat20·Pbat22!+oPbat3=RVocPbat2.Except for SOC, another significant battery state SOE has been accepted extensively [18]. It is more convenient in formulating the theorem and its proof from control design point of view. The battery SOE is defined as the amount of energy stored in battery, divided by the maximum energy capacity of it. The definition of SOE can be written as(6)dSOEdt=-ηbatEmax⁡dEbatdt=-PbatEmax⁡-REmax⁡VocPbat2Emax⁡=Qmax⁡Voc,max⁡,where ηbat is the charge or discharge efficiency of battery and Qmax is the maximum capacity.


## 2.1. Engine Model

The aim of the optimal control strategy of PHEV is to minimize the fuel consumption and deplete the battery capacity over a driving cycle. Fuel consumption of diesel engine is generally modeled as a map for every possible combination of speed and torque. Obviously, the map can only be represented as nonlinear high order function. For simplicity, appropriate Willans line model is usually used to express the function of the engine power and speed. From the characteristic of Willans line model, we can say that the efficiency of the energy conversion device can be modeled by representing the input power as an affine function of the output power and losses.At any given speed, the engine powerPice can be represented as an affine function of the fuel power Pful. The gradient and intercept of each Willans line can be expressed as polynomial functions normally depending on engine speed, by(2)Pfult=e0ωicet+e1ωicetPicet,e0ωicet=e00+e01ωicet+e02ωice2t,e1ωicet=e10+e11ωicet+e12ωice2t,m˙ft=1QLHVe0ωicet+e1ωicetPicet,where e0e1 are the variable coefficient of the fitted Willans model, as shown in Figures 2 and 3. e00,  e01,  e02,  e10,  e11, and e12 are the constants of fitting. The final fitted Willans model is compared with original fuel map in Figure 4.Figure 2
Curve of coefficiente0 with engine speed.Figure 3
Curve of coefficiente1 with engine speed.Figure 4
Comparison of fitted Willans model and original fuel map.


## 2.2. Battery Model

Traditionally, the battery SOC is appropriate for defining the rate change of power energy stored in the battery. The simplified model of battery can be expressed in (3) according to internal resistance model:(3)Ibat=Voc-Voc2-4RPbat2R,where Ibat is the discharge current of battery; Pbat is the power of battery; Voc is the open voltage of battery; R is the internal resistance of battery.The total electric power loss over the battery resistance is(4)PlPbat=RIbat2=Voc-Voc2-4RPbat24R. Assuming the power loss of battery as constant, the power loss can be approximated by formula of expanded Taylor series around Pbat=0, yielding(5)PlPbat≈Pl0+∂Pl∂Pbat0+∂2Pl∂Pbat20·Pbat22!+oPbat3=RVocPbat2.Except for SOC, another significant battery state SOE has been accepted extensively [18]. It is more convenient in formulating the theorem and its proof from control design point of view. The battery SOE is defined as the amount of energy stored in battery, divided by the maximum energy capacity of it. The definition of SOE can be written as(6)dSOEdt=-ηbatEmax⁡dEbatdt=-PbatEmax⁡-REmax⁡VocPbat2Emax⁡=Qmax⁡Voc,max⁡,where ηbat is the charge or discharge efficiency of battery and Qmax is the maximum capacity.


## 3. Pontryagin’s Minimum Principle

It is well known that calculus of variations is restricted to only solve the optimal control problem whose control variable is not constrained. However, the control vector is always under certain constraints in actual engineering system. As a result, classical method of calculus of variations is ineffective in dealing with the problem. On the basis of calculus of variations, PMP proposed by mathematician Pontryagin has become one of the powerful solutions for solving optimal control problem with control vector constrained.In point of view of optimal control, the control strategy of PHEV for specific driving condition can be regarded as the two-point boundary value problem. The initial and terminal time and state are fixed as prior condition. According to PMP, the necessary condition of the optimality for the optimal control is listed below.Without loss of generality, the performance index with constraints of the terminal state is as shown in (7). The state function can be expressed as (8):(7)
J
=
Φ
x
t
f
,
t
f
+
∫
t
0
t
f
L
x
t
,
u
t
,
t
d
t
(8)
s.t.
d
x
d
t
t
=
f
x
t
,
u
t
,
t
.The constraint of the control variable is given by(9)umin⁡<ut<umin⁡,u∈Rn.For PHEV control strategy, (8) represents battery SOE state function expressed by (6). The physical meanings of symbols x,  u,  L, and Φ are battery SOE, power of battery, instant fuel consumption rate of engine, and terminal state requirement expression, respectively.The Hamilton function deduced by (7) and (8) can be formulated as (10)Hxt,ut,λt=Lxt,ut,t+λfxt,ut,t.The state function can be expressed as(11)dxdt=∂H∂λ.The costate can be expressed as(12)dλdt=-∂H∂x.The minimal condition of PMP is(13)Hx∗t,u∗t,λ∗t=minu∈Ru⁡Hx∗t,ut,λ∗t.Initial and terminal state conditions, inner constrains, and transversal condition are formulated as (14)–(16) separately:(14)ϑxt0,t0,xtf,tf=0,ϑ∈Rp,(15)gxt,t≤0,g∈Rp,(16)λtf=∂Φ∂x+∂ϑ∂xtf,where physical expression of g is the constraint of the SOE state variable; that is, SOEmin<SOE<SOEmax.An issue that needs to be paid attention to is that the PMP provides only the necessary condition of optimality instead of sufficient condition. It turned out that several results which meet with the optimality requirements can be worked out with different initial value of costate from the above equations. Previously research has proved that there is only one theatrical globally optimal control strategy for specific driving condition [19]. Decision of the initial value of costate properly becomes the critical step of optimal control of PHEV.


## 4. Costate Estimation Based on Legendre Pseudospectral Method

Proper application of PMP is vulnerable with the initial costate variable. For control strategy of PHEV, the initial costate depends on the driving condition which is difficult to be predicted as described in proceeding. In addition, the optimum of Hamilton function is easy to converge partially for nonlinear singular optimal control problem. Different from PMP which calculate the OCP indirectly, Legendre Pseudospectral Method belongs to the direct OCP calculation method as shooting method. Legendre Pseudospectral Method discretizes the state and control variables at LGL nodes. The variables can be approximated by Lagrange polynomial. Then, the differential equation of state variable and integration equation of cost function can be converted into algebraic operation. Finally, the OCP can be transformed to NLP with control and state variables optimized at LGL nodes. The general framework of costate estimation based on Legendre Pseudospectral Method will be shown as follows.


### 4.1. Time Interval Transformation

The OCP presented in (7)–(10) is formulated over the time interval [t0,tf]. The LGL nodes lie in the interval [-1,1]. For the interval of Legendre polynomial, the following transformation is necessary to express the problem for t∈[t0,tf]⇔τ∈[-1,1]:(17)τ=2t-tf-t0tf-t0,τ∈-1,1.


### 4.2. Legendre Collocation at LGL Nodes

Collocation of Legendre Pseudospectral Method is searching the zero of first-order derivative ofN degree Legendre polynomial L˙N(τ) on the interval. Let τi(i=0,1,…,N) be the LGL nodes. The control and state variables will be discretized in terms of their values at the LGL nodes approximately. Therefore, (N+1) discretized state variables X0,X1,…,XN,Xi∈RN and (N+1) discretized control variables U0,U1,…,UN,Ui∈RN will be achieved. The continuous variables can be approximated by N degree polynomials of the form(18)xτ≈Xτ=∑i=0NLiτXi,uτ≈Uτ=∑i=0NLiτUi,where, for i=0,1,….,N,(19)
L
i
τ
=
∏
j
=
0
,
j
≠
i
N
τ
-
τ
j
τ
i
-
τ
j
,
i
=
0,1
,
2
,
…
,
N
,are the N order Lagrange polynomials. Another characteristic equation of the Lagrange polynomials can be obtained:(20)Liτ=1,i=j0,i≠j.


### 4.3. Transformation of State Equation

Using the discretization processing, the time derivative of the approximated state vector can be approximated by derivative of Lagrange polynomials at the LGL nodes. It can be written as(21)x˙τk≈X˙τ=∑i=0NL˙iτkXi=∑i=0NDkiXi,where k=0,1,…,N; Dki are entries of the (N+1)(N+1)-order differentiation matrix, representing the differentiation of Lagrange polynomials at LGL nodes:(22)
D
k
i
=
L
N
τ
k
L
N
τ
i
τ
k
-
τ
i
,
i
≠
k
;
-
N
N
+
1
4
,
i
=
k
=
0
;
N
N
+
1
4
,
i
=
k
=
N
;
0
,
otherwise
.Similarly, the state equation (10) can be converted to (N+1) constraint equations by collocating at the LGL nodes:(23)∑i=0NDkiXi-t0-tf2fXk,Uk,τk=0.


### 4.4. Transformation of Cost Function

The integral portion of the cost function can be calculated by Gauss-Lobatto integration method. The converted cost function is(24)J=Φxtf,tf+tf-t02∑i=0NωiLXi,Ui,τi,where the weighting factor ωi can be expressed as(25)ωi=∫-11Liτdτ=2NN+1LN2τi.


### 4.5. Converting OCP into NLP

Substituting the discretized equations (18), (23), and (24) into (8)–(10), the OCP can be converted into appropriate NLP, with control and state variables optimized at LGL nodes:(26)J=Φxtf,tf+tf-t02∑i=0NωiLXi,Ui,τis.t.∑i=0NDkiXi-t0-tf2fXk,Uk,τk=0ϑX0,τ0,XN,τN=0,N∈RpgXk,Uk,τk≤0,g∈Rpk,i=0,1,…,N.Further analysis shows that the total number of optimized variables is the sum of discretized control and state number multiplied byN. Excellent solver such as Matlab-fmincon and GPOPS can be utilized for the high dimension sparse matrix of NLP.


### 4.6. Costate Estimation

As introduced in part 1, there is no existing convenient approach of designating the proper costate in real driving condition except for iterative mathematical method. Intuitively first-order necessary optimality conditions can be used to deduce correlation between KKT (Karush-Kuhn-Tucker) multiplier and OCP discretized costate [20]. Specific energy management optimization of PHEV which includes initial and terminal constraints, state equation constraints, will be investigated to justify the correlation.The costate variable equation of Hamilton function can be approximated as(27)λ˙τk=∑i=0NDkiλτi.Imposing the partial derivative of Hamilton function with state variable onto the upper equation, we can deduce(28)tf-t02∂L∂x+∂f∂xλτk+∂g∂xμτk=-∑i=0NDkiλτi.Lagrange function of NLP can be expressed as (29)J~=Φ+tf-t02∑i=0NωiLi+κ~ϑ+∑i=0Nλ~it0-tf2fi-X˙i+μ~igi,where λ~,  κ~, and μ~ are the KKT multiplier of NLP. From the KKT condition, we can deduce the next equation:(30)∂J~∂Xk=0,∂J~∂Uk=0,μ~kgk=0.For the first equation, we can get(31)∂J~∂Xk=tf-t02∂Lk∂Xkωk+∂fk∂Xkλ~k+∂gk∂Xkμ~k-∂∂Xk∑i=0Nλ~iX˙i=0,where(32)∂∂Xk∑i=0Nλ~iX˙i=∑i=0Nλ~i∂∂Xk∑i=0NDinXn=∑i=0NDikλ~i,for (33)Dik=Dki=0,i=k,ωiDik=-ωkDki,i≠k.Therefore, we can find that(34)λτk=λ~kωk,k=1,2,…,N-1.The equation demonstrates the correlation between λ(τk) and λ~k. When k is 0 or N, the effect of Φ and ϑ should be considered to calculate λ(t0) and λ(tN):(35)tf-t02∂L0∂X0ω0+∂f0∂X0λ~0+∂ϑ∂X0κ~+∂g0∂X0μ~0-∑i=0NDi0λ~i=0,as(36)D00=-12ω0,ωiDi0=-ω0D0i,i≠0,∑i=0NDi0λ~i=-ω0∑i=0ND0iλ~iωi-λ~0ω0.Therefore, we can deduce(37)tf-t02∂L0∂X0+∂f0∂X0λ~0ω0+∂g0∂X0μ~0ω0=-∑i=0ND0iλ~iωi-1ω0λ~0ω0+∂ϑ∂X0κ~0. A similar calculation of partial derivative J~ to XN is(38)tf-t02∂LN∂XN+∂fN∂XNλ~NωN+∂gN∂XNμ~NωN=-∑i=0NDNiλ~iωi+1ωNλ~NωN-∂Φ∂XN-∂ϑ∂XNκ~f.For another transverse condition,(39)λ~0ω0+∂ϑ0∂XNκ~0=0,λ~NωN-∂Φ∂XN+∂ϑf∂XNκ~f=0.We can deduce the correlation when k=0 or N:(40)λt0=λ~0ω0,λtN=λ~NωN.As brief summary, the correlation between KKT multiplier of NLP and costate matrix of OCP can be established from (34) and (40). Intuitively, we can use the correlation to estimate the costate of PMP-based control strategy of PHEV.


## 4.1. Time Interval Transformation

The OCP presented in (7)–(10) is formulated over the time interval [t0,tf]. The LGL nodes lie in the interval [-1,1]. For the interval of Legendre polynomial, the following transformation is necessary to express the problem for t∈[t0,tf]⇔τ∈[-1,1]:(17)τ=2t-tf-t0tf-t0,τ∈-1,1.


## 4.2. Legendre Collocation at LGL Nodes

Collocation of Legendre Pseudospectral Method is searching the zero of first-order derivative ofN degree Legendre polynomial L˙N(τ) on the interval. Let τi(i=0,1,…,N) be the LGL nodes. The control and state variables will be discretized in terms of their values at the LGL nodes approximately. Therefore, (N+1) discretized state variables X0,X1,…,XN,Xi∈RN and (N+1) discretized control variables U0,U1,…,UN,Ui∈RN will be achieved. The continuous variables can be approximated by N degree polynomials of the form(18)xτ≈Xτ=∑i=0NLiτXi,uτ≈Uτ=∑i=0NLiτUi,where, for i=0,1,….,N,(19)
L
i
τ
=
∏
j
=
0
,
j
≠
i
N
τ
-
τ
j
τ
i
-
τ
j
,
i
=
0,1
,
2
,
…
,
N
,are the N order Lagrange polynomials. Another characteristic equation of the Lagrange polynomials can be obtained:(20)Liτ=1,i=j0,i≠j.


## 4.3. Transformation of State Equation

Using the discretization processing, the time derivative of the approximated state vector can be approximated by derivative of Lagrange polynomials at the LGL nodes. It can be written as(21)x˙τk≈X˙τ=∑i=0NL˙iτkXi=∑i=0NDkiXi,where k=0,1,…,N; Dki are entries of the (N+1)(N+1)-order differentiation matrix, representing the differentiation of Lagrange polynomials at LGL nodes:(22)
D
k
i
=
L
N
τ
k
L
N
τ
i
τ
k
-
τ
i
,
i
≠
k
;
-
N
N
+
1
4
,
i
=
k
=
0
;
N
N
+
1
4
,
i
=
k
=
N
;
0
,
otherwise
.Similarly, the state equation (10) can be converted to (N+1) constraint equations by collocating at the LGL nodes:(23)∑i=0NDkiXi-t0-tf2fXk,Uk,τk=0.


## 4.4. Transformation of Cost Function

The integral portion of the cost function can be calculated by Gauss-Lobatto integration method. The converted cost function is(24)J=Φxtf,tf+tf-t02∑i=0NωiLXi,Ui,τi,where the weighting factor ωi can be expressed as(25)ωi=∫-11Liτdτ=2NN+1LN2τi.


## 4.5. Converting OCP into NLP

Substituting the discretized equations (18), (23), and (24) into (8)–(10), the OCP can be converted into appropriate NLP, with control and state variables optimized at LGL nodes:(26)J=Φxtf,tf+tf-t02∑i=0NωiLXi,Ui,τis.t.∑i=0NDkiXi-t0-tf2fXk,Uk,τk=0ϑX0,τ0,XN,τN=0,N∈RpgXk,Uk,τk≤0,g∈Rpk,i=0,1,…,N.Further analysis shows that the total number of optimized variables is the sum of discretized control and state number multiplied byN. Excellent solver such as Matlab-fmincon and GPOPS can be utilized for the high dimension sparse matrix of NLP.


## 4.6. Costate Estimation

As introduced in part 1, there is no existing convenient approach of designating the proper costate in real driving condition except for iterative mathematical method. Intuitively first-order necessary optimality conditions can be used to deduce correlation between KKT (Karush-Kuhn-Tucker) multiplier and OCP discretized costate [20]. Specific energy management optimization of PHEV which includes initial and terminal constraints, state equation constraints, will be investigated to justify the correlation.The costate variable equation of Hamilton function can be approximated as(27)λ˙τk=∑i=0NDkiλτi.Imposing the partial derivative of Hamilton function with state variable onto the upper equation, we can deduce(28)tf-t02∂L∂x+∂f∂xλτk+∂g∂xμτk=-∑i=0NDkiλτi.Lagrange function of NLP can be expressed as (29)J~=Φ+tf-t02∑i=0NωiLi+κ~ϑ+∑i=0Nλ~it0-tf2fi-X˙i+μ~igi,where λ~,  κ~, and μ~ are the KKT multiplier of NLP. From the KKT condition, we can deduce the next equation:(30)∂J~∂Xk=0,∂J~∂Uk=0,μ~kgk=0.For the first equation, we can get(31)∂J~∂Xk=tf-t02∂Lk∂Xkωk+∂fk∂Xkλ~k+∂gk∂Xkμ~k-∂∂Xk∑i=0Nλ~iX˙i=0,where(32)∂∂Xk∑i=0Nλ~iX˙i=∑i=0Nλ~i∂∂Xk∑i=0NDinXn=∑i=0NDikλ~i,for (33)Dik=Dki=0,i=k,ωiDik=-ωkDki,i≠k.Therefore, we can find that(34)λτk=λ~kωk,k=1,2,…,N-1.The equation demonstrates the correlation between λ(τk) and λ~k. When k is 0 or N, the effect of Φ and ϑ should be considered to calculate λ(t0) and λ(tN):(35)tf-t02∂L0∂X0ω0+∂f0∂X0λ~0+∂ϑ∂X0κ~+∂g0∂X0μ~0-∑i=0NDi0λ~i=0,as(36)D00=-12ω0,ωiDi0=-ω0D0i,i≠0,∑i=0NDi0λ~i=-ω0∑i=0ND0iλ~iωi-λ~0ω0.Therefore, we can deduce(37)tf-t02∂L0∂X0+∂f0∂X0λ~0ω0+∂g0∂X0μ~0ω0=-∑i=0ND0iλ~iωi-1ω0λ~0ω0+∂ϑ∂X0κ~0. A similar calculation of partial derivative J~ to XN is(38)tf-t02∂LN∂XN+∂fN∂XNλ~NωN+∂gN∂XNμ~NωN=-∑i=0NDNiλ~iωi+1ωNλ~NωN-∂Φ∂XN-∂ϑ∂XNκ~f.For another transverse condition,(39)λ~0ω0+∂ϑ0∂XNκ~0=0,λ~NωN-∂Φ∂XN+∂ϑf∂XNκ~f=0.We can deduce the correlation when k=0 or N:(40)λt0=λ~0ω0,λtN=λ~NωN.As brief summary, the correlation between KKT multiplier of NLP and costate matrix of OCP can be established from (34) and (40). Intuitively, we can use the correlation to estimate the costate of PMP-based control strategy of PHEV.


## 5. Simulation Results and Discussion

To illustrate the effectiveness of method presented in preceding sections, we establish the longitudinal dynamic model in Matlab/Simulink and implement the costate estimation of PMP based on Pseudospectral Method (CE-PMP for abbreviation) in the model. The model has been simulated in user-defined driving cycle which is composed of 3 FHDS. The driving profile is long enough for the PHEV to deplete the battery at the end of the simulation. Initial and terminal state variables SOE are designated as 0.75 and 0.2, respectively. The specification of the vehicle is described in Table1. For distinct comparison, the representative global optimal indirect resolver DP has been used to calculate the theoretic optimal values. Moreover, in order to decide the proper initial costate value, shooting method has been adopted in this paper. The main idea of shooting method is looking for the proper initial value by iteration in order to keep the terminal SOE at the SOEmin when driving cycle finished.Table 1
Main parameters of the PHEV.
Components
Parameters
Value
Vehicle
Curb weight (kg)
1440
Total weight (kg)
1965
Wind resistance
0.32
Rolling resistance
0.0135
Turbo diesel engine
Displacement (mL)
1400
Rated power (kW/r⋅min−1)
62/4500
Max torque (N⋅m/r⋅min−1)
140/(2000–2500)
Transmission
Ratio (ig)
3.58, 2.02, 1.35, 0.98, and 0.81
Differential ratio (i0)
3.947
Battery
Type
NiMH
Rated voltage (V)
288
Rated capacity (Ah)
65
Motor
Type
PMSM
Rated power (kW)
15Figure5 shows the historic plot of vehicle speed with time, respectively. The outstanding accordance can be found in the results between CE-PMP and DP. Figure 6 shows the historic plot of control variable with time, respectively. Figure 7 represents the time historic graph of costate variable. The marks in the plot mean the LGL allocate nodes, where triangle symbol represents numerous circular LGL nodes. Obviously, the accurateness of calculation achieved can be improved by the increase of LGL nodes. Burden of computational resource is more stringent consequently. From the zoom figure we can see that “tiny” jump condition occurred and the CE-PMP can capture the jump if LGL nodes are allocated sufficiently. The jump is “tiny” because the terminal SOE state almost arrives at lower limit exactly at the end of the trip. By comparison, the battery SOE will trigger the lower limit in advance if costate is designated as constant and the violent jump condition occurred. In this case, the inferior optimal rather than global optimal control strategy can be deduced. Proof of this assumption can be found in [6].Figure 5
Vehicle velocity versus time.Figure 6
Control variable versus time.Figure 7
Costate variable versus time.Figure8 shows us the historic plot of SOE state variable. Simulation results for different initial constant costate are plotted in the same figure for easier contrast. Separately, Figures 9 and 10 show the integrated and instant fuel consumption for different costate in comparison to the CE-EMP. From the results, we can find that the method presented in the paper can deduce the closer value of global optimal value than constant initial costate value. With consideration of the heavy burden of computation, the advantage of the method is that it can provide guidance for the initial costate guess and jump condition approximation for PMP real-time application.Figure 8
State variable versus time.Figure 9
Integrated fuel consumption versus time.Figure 10
Instant fuel consumption versus time.Eventually, we listed the results of different number of LGL nodes compared with DP results in Table2. Five cases represent 2, 3, 4, 6, and 8 nodes per second, respectively (the time period of the 3 FUDS driving cycle is 3061 s). Legendre Pseudospectral Method can converge to the approximated minimum of DP results by at least 4 LGL nodes per second. Moreover, computation time rises enormously with increase of number of LGL nodes.Table 2
Computation time and accuracy of different LGL nodes.
Number of nodes
Computation time/s
J
m
i
n/g
error/g
6122
462
2596.2
1.8913
e
3
9183
1072
1572.71
0.867
e
3
12244
2796
770.60
65.75
18366
4932
741.35
36.50
24488
10932
729.46
24.61
DP
1066
704.85
/


## 6. Conclusion

In the paper, the OCP problem of PHEV has been converted to NLP problem with Legendre Pseudospectral Method. From the discretized optimization, estimated costate can be approximated as the Lagrange multipliers of NLP divided by the LGL weights. Simulation results indicate that the costate estimated by our method is closer to those obtained by DP method than constant costate. The paper proposes a kind of general approach of costate estimation for predefined driving condition. Moreover, the method can provide guidance for the initial costate guess and jump condition approximation for PMP application.



---

*Source: 9019450-2016-03-27.xml*