Torque demand control by nonlinear MPC for speed control of vehicles with variable valve lift engine

6th IFAC Symposium Advances in Automotive Control Munich, Germany, July 12-14, 2010

Torque demand control by nonlinear MPC for speed control of vehicles with variable valve lift engine Y. Sakai ∗ M. Kanai ∗ M. Yamakita ∗ ∗

Department of Mechanical and Control Eng., Tokyo Institute of Technology, Tokyo,Japan (Tel: 03-5734-3307; e-mail: [email protected]).

Abstract: In this paper, a nonlinear Model Predictive Control (MPC) is applied for a torque demand control for speed control of vehicles with variable valve lift engine. In the control several constraints, e.g., A/F ratio, valve lift range, etc. are considered. The validity of the proposed method will be demonstrated by numerical simulations using a SICE engine benchmark simulator. Keywords: Engine control; Torque control; Model based control; Nonlinear control. 1. INTRODUCTION Recently, to solve the energy issue and the global environmental problem, auto-makers have developed a high functional engine such as a variable valve lift system. While the development period is becoming shorter to meet needs of the market, developments of engine control systems are also becoming more complex. In this situation, it is difficult to develop products satisfying safety, handle-ability, etc, in the short term. As one of promising directions to resolve this issue, Model Based Development(MBD), which is a development approach using plant models and simulators, attracts much attentions. To support the realization of MBD environments, the SICE benchmark problem was provided (Ohata et al. [2007]). In this study, we apply nonlinear model predictive control(MPC) since the engine system is a MIMO system and requires input constraint consideration. Although MPC has been applied to process controls such as chemical industry conventionally, real time algorithms are developed for fast dynamics system such as mechanical systems (Ohtsuka [2004], Murayama and Yamakita [2009]), recently. In this paper, we apply MPC to an engine system which has a comparatively short sampling period and achieve the control specifications. Since it is becoming more difficult to capture the entire picture of an automobile development due to specialization of the works, each of the devisions requires to use torque as an interface of the development (Gerhardt et al. [1998]). With the background like this, in this paper, we use an approach of torque demand control, which attains car speed tracking control by calculating desired torque from desired car speed, and we use nonlinear MPC for torque tracking control for the desired torque. This paper is organized as follows: Section 2 gives the bench mark model and control specification. In Section 3 we derive the control model of an engine and a car body. Section 4 contains explanations about a procedure to 978-3-902661-72-2/10/$20.00 © 2010 IFAC

494

calculate desired torque, nonlinear MPC, and fuel injection control. Numerical simulations show the effectiveness of these method in Section5. Section 6 is devoted to a conclusion and some ideas of our future work. 2. PROBLEM SETTING A conceptual diagram of an engine-vehicle model is shown in Fig. 1. The plant consists of an engine, an automatic transmission and a car body. We deal with a V6 SI engine with 16-inputs and 5-outputs, which is described as a mixed system of continuous and discrete time system of the following equations: x(t) ˙ = f (x(t), u(t, z)), y(t) = g(x(t), u(t, z)), where x(t) ∈ R , u(t, z) ∈ R16 , y(t) ∈ R5 and z denote state, input, measurable output and engine cycle, respectively. nx denotes the number of the state and should be nx

Fig. 1. The benchmark model Table 1. Notation of input and output uth uvti uvte uvl uf i usa

throttle angle inlet valve phase angle exhaust valve phase angle inlet valve lift fuel amount of cylinder ignition timing of cylinder

yrpm engine speed yst crank angle ymt throttle air mass flow yaf air fuel ratio

10.3182/20100712-3-DE-2013.00164

AAC 2010 Munich, Germany, July 12-14, 2010

determined by modeling the system. The components of input u(t, z) and output y(t) are shown in Table. 1. uth (t), uvti (t), uvte (t) and uvl (t) can be operated continuously, but uf i [z] and usa [z] can be operated only one time in each engine cycle. And the upper and lower bound constraints, 0 ≤ uth ≤ 90[deg] and 0 ≤ uvl ≤ 1, exist due to the limits of the actuators. The output signal yaf cannot be measured in engine starting time because the air fuel ratio sensors do not activate unless the temperature exceeds a certain value. The engine torque τe is determined by the engine speed, the cylinder air mass, the air fuel ratio and the ignition timing, etc. and transmitted to tyres through torque converter, automatic transmission and drive shaft, etc. The engine speed is derived from torque transmitted to tyres, and from it, vehicle speed is determined. The control specifications in the benchmark problem are given as follows: • torque demand control: The engine control system which attains torque tracking control should be constructed after calculating the desired torque based on the desired vehicle speed. • vehicle speed tracking control: The vehicle speed V (t) should follow the desired speed Vr (t). Note that the acceleration should be smaller than a certain value not to aggravate driveability. • fuel consumption: fuel consumption uf which is used to realize vehicle speed tracking should be minimized. 3. MODELING FOR CONTROLLER DESIGN We will now describe a plant model needed to design the control system. We derive a physical model from input to vehicle speed by analyzing the engine simulator firstly, and construct a model for controller design by approximating and discretizing it. In this paper, we focus on deriving a detailed model as possible to obtain an optimized control profiles. 3.1 Engine Model Here, we derive physical models from the input to the engine speed. A conceptual diagram is shown in Fig. 2.

Inlet System Model The dynamics of the air mass inside the intake manifold Mm and the pressure inside the manifold Pm are derived by the conservation law of energy and mass as follows:       6 (M˙ iv )i M˙ th X d Mm κ − 1 , (1) = κ − 1 ˙  − Eth (E˙ iv )i dt Pm i=1 Vm Vm where Vm and κ denote the volume of the intake manifold and the specific heat ratio respectively, M˙ th and E˙ th are the throttle air mass flow and the energy flow, and M˙ iv and E˙ iv are the inlet valve air mass flow and the energy flow. The dynamics of the air mass inside each cylinder Mcyl and the pressure inside each cylinder Pcyl can be derived in a similar way. However, note that switching of the strokes forces the dynamics of Mcyl and Pcyl change.   M˙ iv − M˙ ev  κ − 1   d Mcyl  (q − q− + E˙ iv − E˙ ev  =  Vcyl + (2) , dt Pcyl   dVcyl κ Pcyl − κ−1 dt where Vcyl is the volume of the cylinder, q+ is the energy of combustion and q− is the heat loss calculated from Woschini’s equation. M˙ ev and E˙ ev denote the exhaust valve air mass flow and energy flow. Fuel Model In the fuel model, the amounts of port and valve residual fuel and the amount of the intake fuel into cylinders are determined at crank angle θ = 360[deg] in every cycle and are given by as follows (Guzzella and Onder [2004]): Fwp [z] = Pp Fwp [z − 1] + Rp uf i [z],

(3)

Fwv [z] = Pv Fwv [z − 1] + Rv uf i [z],

(4)

fc [z] = (1 − Pp )Fwp [z − 1] + (1 − Pv )Fwv [z − 1] + (1 − Rp − Rv )uf i [z], (5) where Fwp is the amount of the port residual fuel, Fwv is the amount of the valve residual fuel and fc is the amount of the admitted fuel. Pp , Rp , Pv and Rv are determined by the engine speed ωe (t), the port temperature Tp (t), the valve temperature Tv (t) and the air mass in the previous cycle Mc [z − 1]. The dynamics of Tp and Tv are given by T˙p =

X

(H1 q− − H2 (Tp − Ta )),

(6)

T˙v = H3 q− − H4 (Tv − Ta ), (7) where H1 , H2 , H3 and H4 denote the coefficient of heat transfer and Ta denotes the atmosphere temperature. The intake air mass into each cylinder Mc [z] are described as following equation by integrating the inlet valve air mass flow M˙ iv during one intake stroke: Z 580 Mc [z] = M˙ iv dθ. (8) 360

Combustion Model In the engine simulator, the combustion energy q+ used in the dynamics of air pressure inside the cylinder Pcyl (2) is modeled as

Fig. 2. Engine model 495

AAC 2010 Munich, Germany, July 12-14, 2010

q+ (z + 1, t) = W (usa [z], θ) · MAPq (fc [z], Mc [z]) · ωe , (9) where W (·) is a Wiebe function and MAPq (·) is given as a map. The input usa is determined at θ = 580[deg] in every cycle. Note that we have equation (9) if 9 ≤ Mc [z]/fc[z] ≤ 22. Otherwise, the right hand side of the equation become 0 because of engine miss firing. Crank-shaft Rotation Model From the translatory movement of each cylinder piston and the rotational movement of the crank, the dynamics of the engine rotational speed ωe is described as J(θ)ω˙ e = τe − τf − τtr , (10)   6 X d(Vcyl )i d2 (Vcyl )i 2 τe = [(Pcyl )i − Pa ] ω − mp e , dθ d θ2 i=1 where J(θ) is the inertia moment of the engine, τe is the engine torque, τf is the engine friction torque, τtr is the load torque from the torque converter, Pa is the atmosphere pressure and mp is the mass of the pistons. 3.2 Torque Transfer System and Vehicle Model

In the engine simulator, an automatic transmission is modeled and the torque transmitted from the torque converter to the engine τtr and to the shaft τs are given as maps as follows:

By simplification and discretization explained above, we obtain a model for controller design as xk+1 = xk + f (xk , uk )∆θ, where k denotes the number of step, and we define

T T xk = [Mm , Pm , ωe2 , Mcyl , Pcyl , TvT , TpT , Toil , Liv , v, T T T Fwv , Fwp ] ∈ R34 .

To check the accuracy of the derived model, we examined the behavior of the model in open-loop, without measurable output. The result is shown in Fig. 3. Since we can see the engine speed ωe is estimated with high accuracy, the derived model can be expected to be accurate enough to predict the behavior of the benchmark model.

(11)

4. CONTROL METHOD

(12)

where ωs denotes the shaft rotational speed. By ignoring skidding of the tyre and effect of the torsional stiffness in the drive-train, the shaft rotational speed ωs and the tyre drive force Ft can be described as G G ωs = v, Ft = τts , (13) rt rt where rt is the tyre radius and G is the gear ratio of the differential gear. Here, we describe a motion equation of the vehicle by modeling only the air resistance as Meq v˙ = Ft − Cv 2 , (14) where Meq denotes the equivalent mass including rotating part inertia and C denotes the air resistance.

In this section, we design the torque demand controller for speed control of vehicles. Fig. 4 shows the structure of the controller. The controller consists of the open-loop state estimator, the torque reference generator and two 3000 actual estimated 2500 Engine Speed [rpm]

τs = MAPC (ωs /ωe )MAPT (ωs /ωe )ωe2 ,

Since the plant is a 4-stroke engine system which has 6 cylinders, each of crank angles is different from one next to it by 120[deg]. Therefore the interval to determine the discrete input usa and uf i is also 120[deg]. Because of this, we execute discretization by converting the physical model from time domain to crank angle domain and by considering 120[deg] of crank angle as 1 step of discretization basically. However, since the inlet system has very fast dynamics, each one step is divided into more steps for Mm , Pm , Mcyl and Pcyl .

uk = [uth , uvl , usa , uf i ]T ∈ R4 ,

We derived the physical model from the input to the engine speed in the previous section. Here, we will derive a physical model from the engine speed to the vehicle speed whose desired value will be given as a target trajectory.

τtr = MAPC (ωs /ωe )ωe2 ,

• Approximation of maps: Map functions in the equations are approximated as polynomials by least-square fitting. • Selection of redundant control input: We use input uth , uvl and usa for torque control because the valve timing inputs uvti and uvte seems to have no great influence on torque production as compared with uth , uvl and usa .

2000 1500 1000 500 0 0

5

10

Fig. 3. Estimation of benchmark model behavior

3.3 Derivation of State-Space Expression We need to express the model for controller as a discrete time system for the MPC control approach considered in Section 4. For this reason, we will derive a model in statespace model by simplifying and discretizing the physical model derived in the previous section. Firstly, we introduce following simplification for control calculation. 496

15

Time [s]

Fig. 4. Structure of the controller

AAC 2010 Munich, Germany, July 12-14, 2010

subcontollers. The nonlinear MPC with input constraints is applied for the torque tracking control, and the fuel injection controller is independently designed in order to avoid misfire.

Firstly, we consider the case where the vehicle speed is tracked for desired trajectory, that is, v = vd , v˙ = v˙ d . In this case, desired values can be calculated by using (12)(14) as τsd = a1 v˙ d + a2 vd2 , ωsd = a3 vd , τsd 1 = MAPC (ωsd /ωed ) · MAPT (ωsd /ωed ) · 2 ωsd (ωsd /ωed )2 := MAPCT (ωsd /ωed ), (15) where a1 = Meq rt /G, a2 = Crt /G and c3 = G/rt . (16)

where MAPCT −1 : τs /ωs2 → ωs /ωe is inverse function of MAPCT . By substituting ωed and ωsd for (11), we can calculate the desired load torque from the torque converter τtrd as 2 τtrd = MAPC (ωsd /ωed ) ωed . (17) Although τtrd is not the engine output torque, we can design the tracking controller for τtrd only with knowledge of the the torque converter property such as (11). Thus we consider the control problem for tracking τtrd as a benchmark problem. However, equation (17) is considered only about the steady state, and the transient torque tracking error will remain as a steady state offset vehicle speed error. In order to reduce the offset error, we compensate τsd using vehicle speed error e = v − vd as    Z a2 2 τ˜sd = a1 v˙ d + v − kp e + ki e dt . (18) a1 By assuming perfect torque tracking is achieved, that is τs = τ˜sd , the vehicle speed error dynamics will be given as Z 0 = e˙ + kp e + ki e dt, (19) where kp , ki is the positive constants. From (19), vehicle speed error e → 0 if τs = τ˜sd , and we can calculate torque reference with error compensation τ˜trd by substituting τ˜sd for (16) and (17). 4.2 Nonlinear MPC We will now describe the nonlinear MPC method for a torque tracking control. The nonlinear MPC problem we address here has following form: N −1 X min ϕ[xN (t)] + L[xk (t), uk (t)]∆τ, (20) uk

λk (t) = λk+1 (t) + HxT [xk (t), λk+1 (t), uk (t)]∆τ, λN (t) = ϕTx [xN (t)], Hu [xk (t), λk+1 (t), uk (t)] = 0.

4.1 Torque Reference Generator

From (15), ωed is represented by ωsd , ωed = 2 ) MAPCT −1 (τsd /ωsd

Let λ denote the costate defined by λ = ∂J/∂x and the Hamiltonian defined by H = L + λT f . The necessary condition of the optimal control is obtained as follows:

k=0

subject to : xk+1 (t) = xk (t) + f (xk (t), uk (t))∆τ, (21) where the sampling time is ∆τ = T /N , xk (t) stands for the predictive value of k step and xd stands for the target value. We calculate the input sequence every step U (t) = {u0 (t), u1 (t), · · · , uN −1 (t)}. 497

N

N

For a given x0 (t), U (t), {xi (t)}i=0 and {λi (t)}i=0 are calculated. Then the necessary condition is regarded as an simultaneous equation as   HuT [x0 (t), λ1 (t), u0 (t)]   HuT [x1 (t), λ2 (t), u1 (t)]   F [U (t), x(t), t] =   = 0. ..   . HuT [xN −1 (t), λN (t), uN −1 (t)]

N The solution of this equation {Ui (t)}i=0 minimizes the function(20). The input value u0 (t) which is only the N step of {Ui (t)}i=0 is applied to the plant.

(22) cost first

Nonlinear MPC servo control Conventionally, this problem can be solved in real time by C/GMRES method proposed in Ohtsuka [2004]. A large linear equation is derived by the continuation method and the linear equation is solved efficiently by GMRES method. However, this method can be applied for regulator problem, and it is needed to calculate input value for steady state(Aoki et al. [2008]). For simple system, the input value of steady state can be calculated by inverse model but for complex models, it is difficult to obtain the input value in advance. Therefore we apply the nonlinear MPC method proposed in Murayama and Yamakita [2009]. In this method, the integrator is introduced in front of the system because a servo control system is required for tracking control. This method is often used in process control, but here, we also consider improvement of computational speed. By the modification, new input v = u˙ is calculated and input to the engine u is treated as state variables. In the performance function, weights are imposed to not u but v. Under such conditions, the necessary condition of the optimal control is calculate as   λ1(nx +1) + R(1,1) v0(1) λ1(nx +2) + R(2,2) v0(2)    = 0, F = ..   . λN (nx +nu ) + R(nu ,nu ) vN (nu ) where λk(n) represents the n-th element of λk and R(i,j) represents the element of R in the i-th row and the jth column. This equations are linear equations with one unknown variable so it is easy to solve input sequence analytically and we don’t have to use iteration method such as Newton method. Furthermore, in the engine system, there are several constraints, e.g., throttle angle, valve lift range, A/F ratio, etc. Here, we consider about input constraints such as umin ≤ u ≤ umax . However, in the servo control configuration, input u is treated as state variable so it is necessary to formulate as state constraints which are hard to be deal with. To avoid this problem we introduce a sigmoid function given by

AAC 2010 Munich, Germany, July 12-14, 2010

5. SIMULATION

Fig. 5. The nonlinear MPC servo control with input constraints

In this section, we apply the proposed method to the SICE engine benchmark problem. Here we consider the torque demand control for speed control of vehicles of benchmark simulator. A target vehicle speed is a first step of 10-15 driving mode. Here, the future target speed is not given.

The sigmoid function is introduced to the system as shown in Fig. 5. In doing so, input variable to the objective system satisfies the constraint condition automatically. Let xu denote the input to the real process, state equation about xu is shown by

The engine benchmark model has a variable valve lift system so engine torque is related to throttle angle, inlet valve lift, inlet and exhaust valve timing, etc. These redundant inputs should be effectively used for improvement of response, reduction of fuel consumption, etc. We use throttle angle, inlet valve lift and ignition timing as input for the torque demand control. The ignition timing input is used for preventing overshoot of the vehicle speed. Pumping loss is expected to be reduced by the combined use of throttle and inlet valve.

dxu ∂xu ∂u = dt ∂u ∂t a (24) = (xu − β)(α + β − xu )v. α Because this state equation is primary expressed by v, the necessary condition of the optimal control can be solved analytically as the case without constraints described above.

In order to apply the nonlinear MPC method, we define the cost function in (20) as 1 ϕ = (yN − yd )T Qf (yN − yd ), (27) 2 1 1 (28) L = (yk − yd )T Q(yk − yd ) + v T Rv, 2 2 where yk is the output variables for evaluation, and yd is the desired trajectory.

However, in this method, there is the problem that sensitivity of the new input v for the cost function becomes too small when input to sigmoid function u becomes so large. In such case, input to the real process xu doesn’t change from the neighborhood of constraints. We can avoid this problem by introducing weight function as   1 1 Qb (xu ) = a + −b . (25) umax − xu xu − umin

In order to achieve torque demand control while considering also the driveability and the fuel consumption, we choose four output variables:

ςa (x) =

1 . 1 + exp(−ax)

(23)

The weight function shown in (25) prevents input xu stacked near the bounds because it will be quickly large when xu approach to the neighborhood of the constraints.

4.3 Fuel Injection Controller The A/F ratio in the cylinder is defined as α = Mc /fc where Mc is air mass and injected fuel fc , and the ideal value is αd := 14.5. With the ideal A/F ratio, not only efficient engine combustion is realized but also three-way catalyst converters work best to avoid the un-desirable exhaust. Therefore in this paper, we calculate the fuel injection input to make the ideal A/F ratio using inverse model of (3)-(5). ˆ c from We can estimate the air mass into the cylinder M (8), and then, we have the fuel injection input to make the ideal A/F ratio as fcd [k] − (1 − Pp )Fwp [k − 6] − (1 − Pv )Fwv [k − 6] uf i [k] = , 1 − Rp − Rv (26) ˆ c /αd is desired where k stands for sampling index, fcd := M injected fuel mass. Although Pp ,Rp ,Pv , and Rp depends on ωe ,Tp ,Tv , and Mc , we can estimate these state variables by the engine model derived in Section 3. 498

• the engine speed ωe : When the gear is ’Neutral’, we consider the engine speed control. To avoid the overshoot of the engine speed, the engine speed is regulated to desired value ωed by the nonlinear MPC. • the load torque from the torque converter τtr : The torque demand control for the load torque τtr is considered when the gear is ’Drive’. The desired value is calculated by the torque reference generator, and the torque tracking is achieved by the nonlinear MPC. • the intake air mass into the cylinder Mc : Here, the fuel injection input is calculated by (26). So, we can reduce the fuel consumption by regulating the intake air mass Mc . • the spark angle xusa : To achieve high efficiency torque generation, we evaluate the difference between the spark angle input and an optimal (Minimum spark advance for Best torque : MBT) spark angle. According to the definition above, output and its desired are defined as y := [ ωe , τtr , Mcyl , xusa ]T , yd := [ ωed , τtrd , 0, φMBT ]T .

(29)

The input to the real process are bounded as T

T

umin := [ 5◦ , 0.01, −30◦ ] , umax := [ 90◦ , 1.0, 30◦ ] . Simulation results are shown in Fig. 6 - Fig. 8. In Fig. 6, we plot the result of engine speed control. The engine speed is regulated to ωed = 1000[rpm] while the

AAC 2010 Munich, Germany, July 12-14, 2010

2

1

0

500

1 Fuel injection [g]

1000

1.2

25 Air−Fuel ratio #1 [−]

1500

x 10

30

desired actual

3

Acceleration [km/h/s]

2000 Engine Speed [rpm]

−4

4

2500

20 15 10 5

0 0

5

10

−1 0

15

5

10

15

0 0

Time [s]

Time [s]

(a) Engine speed

0.8 0.6 0.4 0.2

5

10

15

0 0

5

Time [s]

(b) Vehicle acceleration

(a) Air-Fuel ratio

Fig. 6. Simulation result: Engine speed control

10

15

Time [s]

(b) Fuel injection

Fig. 8. Simulation result: A/F control 6. CONCLUSION

100

desired actual

desired actual

20 Velocity [km/h]

Torque from TM to Engine [N.m]

120

80 60

15

10

40

5 20 0 0

5

10

0 0

15

5

10

(a) Load torque from torque-converter to engine

(b) Vehicle speed

Future work will be focused on both more control oriented modelling and developments of MPC. Especially, to avoid the misfire at engine startup, we consider applying the nonlinear MPC to be able to handle state constraints such as the collocation method in Pytlak [1999].

1.2

100

1 Valve Lift [−]

80

60

40

20

0 0

15

Time [s]

Time [s]

Throttole angle [deg]

In this paper, torque demand control method by nonlinear MPC has been proposed. To deal with torque demand control problem, we design the controller which consists of the torque reference generator and nonlinear MPC. The validity of the proposed method is demonstrated by numerical simulations using a SICE engine benchmark simulator. From the simulation results, it has been shown that the proposed method satisfy the design specification while input constraints are satisfied.

25

140

0.8 0.6 0.4 0.2

5

10

0 0

15

5

10

REFERENCES

15

Time [s]

Time [s]

(c) Throttle angle

(d) Valve lift

Ignition timing [deg]

50

0

−50 0

5

10

15

Time [s]

(e) Ignition timing

Fig. 7. Simulation result: Torque demand control

gear is ’Neutral’(t ≤ 3[sec]), and it prevents the drastic change of vehicle acceleration. Fig. 7 shows the result of torque control. At first, from Fig. 7(a), we can see that torque tracking is achieved by nonlinear MPC, and torque reference is modified when vehicle speed error is occurred. As a result, the vehicle speed follows a target speed closely. Furthermore, in Fig. 7(c)(e), it also be seen that input constraints are satisfied. The A/F control result are shown in Fig. 8. From Fig. 8(a), we can see the A/F ratio stays near the ideal A/F in the steady state. Fig. 8(b) shows that the fuel injection input. The minimization of fuel consumption is considered indirectly by minimizing the intake air mass into the cylinder and the difference between the spark angle input and the MBT spark angle. 499

J. Aoki, T. Ohtani, and T. Ohtsuka. A Compensation Method of Offset Caused by Nonlinear Model Predictive Control using Continuation/GMRES Method (in Japanese). In Proceedings of Joint Conference on Automatic Control, pages 910–915, 2008. J. Gerhardt, H. H¨onninger, and H. Bischof. A New Approach to Functional and Software Structure for Engine Management Systems - BOSCH ME7. SAE Paper, (980901), 1998. L. Guzzella and C. H. Onder. Introduction to Modeling and Control of Internal Combustion Engine Systems. Springer-Verlag, 2004. A. Murayama and M. Yamakita. Speed Control of Vehicles with Variable Valve Lift Engine by Nonlinear MPC. In Proceedings of SICE 2009, 2009. A. Ohata, J. Kako, T. L. Shen, and K. Ito. Benchmark Problem for Automotive Engine Control. SICE Annual Conference, 2007. T. Ohtsuka. A continuation/GMRES method for fast computation of nonlinear receding horizon control. Automatica, 40:563–574, 2004. R. Pytlak. Numerical Methods for Optimal Control Problems with State Constraints. Lecture Notes in Mathematics. Springer Verlag, 1999.