Modeling and control of a parallel waste heat recovery system for Euro-VI heavy-duty diesel engines

.


Introduction
Nowadays, the primary power source for transportation is provided by internal combustion engines.Despite efforts to improve fuel economy in modern engines, approximately 60%-70% of the fuel power is still lost through the coolant and exhaust [1].Moreover, the upcoming CO 2 emission legislation needs to be fulfilled.For heavy-duty vehicles, USA legislation indicates a 20% CO 2 emissions reduction by 2020.In Europe, similar requirements are expected to be introduced.To meet these upcoming CO 2 regulations and further reduce fuel consumption, diesel engines equipped with waste heat recovery (WHR) systems based on an organic Rankine cycle (ORC) seem very promising [2][3][4], especially for long-haul heavy-duty truck applications [5,6].A WHR system allows energy to be recovered from the exhaust gas to generate mechanical power useful for propulsion.
For automotive applications, control of WHR systems is challenging.This is due to highly dynamic engine conditions, to interaction between the engine and the WHR system and to constraints, namely actuator limitations and safe operation.The control strategies dedicated to automotive applications are mostly based on a proportional-integral-derivative (PID) control type of approach [7][8][9][10][11][12].However, this approach is most well established for single-input single-output systems.For multiple-input multiple-output systems with constraints, a model predictive control (MPC) can be developed, which, to the authors knowledge, has not been used for automotive WHR applications.Predictive control has been applied only to power plants [13,14] and refrigeration systems [15].The main difference is that power plants and refrigeration systems work at different time scales as compared to automotive WHR systems that are governed by highly dynamic engine behavior.
In this paper, a switching model predictive control strategy is presented that considers the effect of highly dynamic engine disturbances caused by real on-road driving conditions.The WHR system is highly non-linear, mainly due to the two-phase flow phenomena.Thus, a global controller can give limited performance over the complete engine operating area.To account for this, a switching mechanism is developed that assigns a controller for a certain engine operating area.As a result, the control strategy allows operation close to the safety limit with good disturbance rejection properties.To benchmark the MPC strategy performance, a classical proportional-integral (PI) control strategy [16] is used, for comparison.Furthermore, an offline nonlinear MPC (NMPC) strategy is implemented to compare the MPC solution with the optimum solution in the case that the WHR system model is exactly known.The NMPC scheme is developed only for comparison purposes, since it is characterized by high computational complexity, making it unsuitable for embedded applications using off-the-shelf hardware.
For control design and optimization, dynamical modeling of WHR systems plays an important role.In [17], a dynamical model of the complete WHR system has been developed.The model captures the two-phase flow phenomena based on the conservation principles of mass and energy.However, in [17], simple evaporator models were used with simplified heat transfer coefficients.To improve the accuracy of the WHR system model, more advanced evaporator models are included, based on the approach used in [18] for the exhaust gas recirculation (EGR) evaporator.This approach combines the finite difference modeling approach with a moving boundary one.Thus, it is able to capture multiple phase transitions along a single pipe flow.
The paper is organized as follows.In Section 2, the experimental set-up is described.In Section 3, the improved WHR system model is presented and validated against the measurements.In Section 4, the control design and switching mechanism are presented.The control performance is evaluated in Section 5 by means of simulation results and comparison with a classical proportional-integral (PI) and NMPC control strategy.Conclusions and directions for future research are discussed in Section 6.

Experimental Set-Up
The experimental set-up is schematically illustrated in Figure 1.A six-cylinder Euro-VI European emission standards heavy-duty diesel engine is studied with an aftertreatment system that consists of: a diesel oxidation catalyst (DOC), diesel particulate filter (DPF), selective catalytic reduction (SCR) and an ammonia oxidation catalyst (AMOX).The engine is furthermore equipped with a WHR system that recovers heat from both the EGR line and main exhaust (EXH) line.The studied WHR system is based on a typical organic Rankine cycle (ORC) with pure ethanol as the working fluid.Ethanol is selected because of its physical and thermodynamic properties, which are suitable for this low temperature application.The ethanol is pumped from a reservoir at atmospheric pressure through plate-fin counter-flow evaporators, where it changes its state from liquid to two-phase and then vapor at the evaporator outlet.By expanding the vaporized ethanol, mechanical power is generated at the expander shaft.In this set-up, a two-piston expander is used that, according to the manufacturer data, can withstand ethanol in the two-phase state for short periods of time.To close the cycle, a condenser is used that brings the ethanol back to the liquid state.For safety reasons, a pressure relief valve is used that limits the system pressure to 60 bar.The measurement system consists of thermocouples to measure temperature and thin-film pressure sensors.The exhaust gas mass flow rate is measured using Coriolis flow meters, while the ethanol mass flow rate is measured using helicoidal flow meters.To test the system under various engine loads and rotational speeds, a dynamometer is connected to the engine.To minimize heat losses to the environment, the WHR system was insulated.
Note that in the studied system, the expander and the pumps are mechanically coupled to the engine crankshaft.A disadvantage of this design is that due to mechanical coupling, the pressure in the system can be restricted, and thus, the power of the WHR system becomes limited.The advantage is that the mechanical power is directly transmitted to the engine, without energy losses from a power converter or variable transmission.During gear shifting, when the requested power is less than the WHR system net output power, P req < P whr , a throttle valve is closed (u t = 0%) to avoid unwanted torque response.However, in this study, we focus only on the power mode (excluding gear shifting), and thus, the throttle valve is kept fully opened (u t = 100%).
To avoid droplets, i.e., to guarantee the vapor state after the EGR evaporator and EXH evaporator, two bypass valves, u egr and u exh , are used.These bypass valves manipulate the ethanol flow rate, such that the vapor state is maintained after the evaporators in the presence of engine disturbances.On the exhaust gas side, the EGR valve u g 1 is controlled using the standard calibration for a conventional engine.From the cooling perspective, a valve u g 2 is used to bypass the exhaust gases at the EXH evaporator and, so, accommodates the maximum condenser cooling capacity.

Waste Heat Recovery System Modeling
The WHR system model is described using a component-based approach using first principle modeling.It consists of a reservoir, pumps, valves, evaporators, condenser and an expander model.The pumps and expander are map-based components, while the rest is modeled based on the conservation principles for mass and energy.The pipes, on the vapor side (indicated with red in Figure 1) are modeled as a compressible volume.The WHR system model is developed based on the following assumptions: • the working fluid flows always in the positive direction; • transport delays and pressure drops along the pipes are neglected; • pressure drops along the evaporators are neglected on both the working fluid and exhaust gas side; • the exhaust gas change in density along the evaporators is neglected; • the condenser is ideal, so that the reservoir provides working fluid at the ambient pressure of 1 bar and a temperature of 65 • C.
Here, a brief review is given for the main components within the system, and more details can be found in [17,18].

Pump
The WHR system is equipped with two identical pumps.The mass flow rate through a pump is computed based on the ideal mass flow rate and volumetric efficiency as: where the volumetric efficiency η vol = f (ω p , p r ) is a function of pump speed ω p (rev/s) and the pressure ratio across the pump p r , according to the measurements.The ideal mass flow rate is given by: where V d is the pump displacement volume and ρ is the working fluid density.The pump power is obtained using: where η is is the isentropic efficiency and accounts for internal pump energy loss caused by friction and external thermal loss.Due to the low viscosity of ethanol, the isentropic efficiency is assumed at the upper limit [19], i.e., η is = 0.9.The outgoing working fluid temperature is calculated based on the energy balance equation: where c p = f (T p,in , p in ) is the specific heat capacity of the working fluid.

Valve
The valves within the WHR system are modeled as incompressible flow valves for flows in the liquid state and compressible flow valves for the two-phase and vapor flows.The thermodynamic state of the outgoing fluid is calculated assuming an isenthalpic process through the valve, i.e., h out = h in .

Incompressible Flow Valve
The mass flow rate through the valve for incompressible flow is expressed as: where u v represents the duty cycle from 0 to 1 applied to the valve actuator, c d is the discharge coefficient, A 0 is the orifice area and p in − p out is the pressure difference across the valve.

Compressible Flow Valve
The compressible flow valve is modeled using: where φ is the compressibility coefficient defined as: , with: The parameter ψ has the following expression, whether the flow is subsonic or supersonic: where p out /p in is the pressure ratio across the valve and 2 γ+1 γ γ−1 is the critical pressure ratio.

Evaporator and Condenser
The thermal behavior of the WHR system is mainly influenced by the evaporators and condenser.Therefore, the modeling of these components plays a crucial role to describe the WHR system dynamic behavior [11].The evaporators are plate-fin counter-flow heat exchangers in which the working fluid changes its state from liquid to two-phase and from two-phase to vapor [8].Compared to [17], an improved model for the EGR evaporator is described in [18].Here, we incorporate the enhanced EGR evaporator model and derive an EXH evaporator using a similar representation as in [18].The enhanced evaporator model combines a finite difference modeling approach with a moving boundary one.Thus, it captures multiple phase transitions along a single pipe flow.To characterize the single-phase flow, as well as the two-phase flow phenomena, more advanced heat transfer coefficients are used, which improve the prediction accuracy of the exhaust gas and working fluid temperature in general.The evaporator model is based on the conservation principles of mass and energy.The time derivative of pressure is neglected, since the pressure dynamics are much faster than the thermal phenomena.As a result, the evaporator model can be written as the following set of partial differential equations: Conservation of mass (working fluid): Conservation of energy: Conservation of energy at the wall: where z is the space coordinate along the length L of the evaporator and V f , V g and V w are the volume occupied by the working fluid, by the exhaust gas and by the wall, respectively.Furthermore, S is the surface area, c p is the specific heat capacity, λ w represents the wall thermal conductivity and δ w is the wall thickness.A constant factor γ w = 0.6 is introduced in Equation ( 12) to separate the mass of the wall, i.e., M w = ρ w V w , including fins on the exhaust side from the mass of the wall on the working fluid side.The wall conduction has been found to be of great importance for the dynamic evaporator behavior.Thus, in Equation ( 12), the transversal wall conduction is computed based on two temperature points, T wg and T wf .
The heat transfer coefficients on the exhaust gas side α g and working fluid side α f are based on Nusselt numbers correlation selected from [20].The heat transfer on the exhaust gas side is enhanced by fins, and thus, a Nusselt number correlation for finned surfaces has been selected.On the working fluid side, the heat transfer coefficient depends on the flow regime, laminar or turbulent, indicated by the Reynolds number, as well as the working fluid state of matter, i.e., liquid, two-phase and vapor.In the two-phase state, the heat transfer coefficient is a function of the vapor fraction obtained from the temperature-enthalpy diagram illustrated in Figure 2. The working fluid vapor fraction is computed as: where h l (p f ) and h v (p f ) denote the specific saturated liquid and vapor enthalpy, respectively, as a function of system pressure p f .Using Equation ( 13), the working fluid has the following states: The condenser model is based on a simplified approach described in [17].

Mixing Junction
The mixing junction is considered to be characterized by fast fluid dynamics, such that the stationary conservation laws for mass flow and energy can be applied.Thus, the outlet mass flow is given by: where n is the total number of inlet fluid streams and ṁk,in denotes the k-th inlet fluid mass flow.Assuming ideal mixing, the stationary mixing junction vapor fraction is calculated with the energy balance equation as: where χ k,in are a function of specific enthalpy and system pressure according to Equation (13).

Pressure Volume
Let us consider the piping after the evaporators as a fixed volume.A pressure volume model is introduced to characterize the pressure and temperature dynamics inside the volume, as a function of the inlet and outlet conditions.The pressure volume model assumes that the working fluid is in the vapor state and that it behaves like an ideal gas.The laws for the mass conservation Equation (16a) and energy conservation Equation (16b) are then applied to this ideal gas: where m is the working fluid mass and T is the temperature inside the pressure volume.
Next, assume the ideal gas law: and define: Then, Equation (16a) and (16b) become: where V is the volume, R is the ideal gas constant and γ is given by Equation (8).For steady-state conditions, Equation (19a) and (19b) simplify to ṁout = ṁin and T = −(b/a)T in = T in , respectively.

Expander
The expander model is a map-based component, in which the outgoing mass flow rate is: where ω exp is the expander speed and p in /p out is the pressure ratio.In a similar way, the WHR system net output power, including the pumps' power, is given by: Next, the expander power is calculated by adding to Equation ( 21) the absolute net power of the pumps: where P p1 , P p2 are the pumps' requested power.Assuming that the heat losses to the environment are negligible, the enthalpy at the expander outlet is obtained from: Similar to the vapor fraction modeling for the evaporator, the expander outlet vapor fraction is obtained by substituting Equation (23) into Equation (13).

Experimental Validation
The WHR system model validation is performed over a wide range of engine operating conditions: engine torque between 500-2600 N•m and engine speed between 1200-1850 rpm.The experimental data was first processed to check if the energy conservation principles are satisfied.To this end, data reconciliation techniques have been applied to correct for the mass flow rate measurement inaccuracies (more details are given in [17]).
The measured signals used as input for the WHR system model are: the expander speed, bypass valves duty cycle, exhaust gas mass flow rate and temperature for both the EGR and EXH evaporator.Based on these input signals, the WHR system model predicts the heat flow rate through the system, the working fluid mass flow rate, the system pressure, the temperature and the power.For brevity, here, we show a few results that are related to the improved evaporator models.
In Figure 3, the heat flow through the system, mass flow rate and pressure are compared with the measurements.All of the heat flow rates are expressed relative to the maximum condenser cooling capacity.The working fluid mass flow rates ṁf,egr and ṁf,exh and system pressure p f are expressed relative to the maximum values that the system can reach.To calculate the exhaust gas and working fluid heat flow rate for both evaporators, the following expressions are used: where c pg is the exhaust gas heat capacity.For brevity, in Figure 3a, we show the total exhaust heat flow rate, calculated as: Qg,tot = Qg,egr + Qg,exh + Qloss,egr + Qloss,exh where Qloss,egr and Qloss,exh are the heat losses to the environment (see Equation (11a)).The exhaust gas heat flow rates Qg,egr and Qg,exh are transferred to the working fluid in the form of Qf,egr and Qf,exh , respectively.From the simulation results, we can see that heat flow rates through the system are well captured by the model on both the exhaust gas side and working fluid side, with an error of 2%.Furthermore, the mass flow rates through both evaporators show good behavior.The average error is 3%, as well as the pressure dynamics, with an average error of 2%.The oscillations observed between 2000-4000 s are a consequence of the WHR control algorithm, which manipulates the bypass exhaust gas valve u g 1 due to cooling requirements.As can be seen, after 4000 s, a peak in the total exhaust gas heat flow rate is present that exceeds the condenser cooling capacity.However, this is for a short period of time, and thus, not visible in the working fluid heat flow rate, which could eventually lead to cooling issues.
In Figure 4, the working fluid temperature and vapor fraction after the evaporators and mixing junction are compared with experimental data.The working fluid temperature after the EXH evaporator shows that there are differences as compared to the measurements.This follows from the applied data reconciliation techniques, which are not perfect and still result in some inaccuracies with respect to the energy conservation principle.To overcome this behavior, more accurate mass flow sensors are necessary.Unfortunately, this equipment was not available at the experimental set-up.Furthermore, the instrumentation to measure the vapor fraction was not available.However, in the liquid region (χ ≤ 0) and vapor region (χ ≥ 1), the vapor fraction can be calculated based on the measured temperature and pressure.The two-phase behavior can be validated by observing at which moment the working fluid switches from the vapor state to the two-phase state and back at the evaporator outlet.This behavior can be seen for the EXH evaporator between 5000-7000 s, where the vapor fraction χ f,exh drops significantly due to the low heat input (see the corresponding Qf,exh in Figure 3a).(b) The model prediction error is calculated using the following expression: where y sim is the output to be compared, while y max and y min are the outputs at maximum and minimum engine thermal load, respectively.The resulting average and maximum error for the heat flow rate, mass flow rate, system pressure and temperature are summarized in Table 1.
Table 1.Summary of the model prediction error and improvement compared with data from [17].
Compared to the previous model, we first see that the working fluid temperature prediction error decreased on average by 2%, while the maximum error decreased by 7%.This indicates that both the steady-state and dynamic behaviors of the evaporators are improved.For the mass flow rate, there is no improvement, since the pump and the valve models are kept the same.Furthermore, there is insignificant improvement for the average heat flow rate and system pressure error.However, due to the improved EGR and EXH evaporator models, the maximum heat flow rate and system pressure error decreased by 2% and 1%, respectively.

Control Design
In this section, the control objective is presented, followed by the model for control and the problem formulation in the model predictive control framework.

Control Objective
Based on Equation ( 13), let us denote the vapor fraction after the EGR and EXH evaporator with χ f,egr and χ f,exh .Using Equation ( 15), the vapor fraction downstream of the mixing junction is given by: In Figure 5, the control scheme is illustrated.The control input is u = [u egr u exh ] ; the output is y = [χ f,egr χ f,exh ] ; and the engine disturbance w = [ Qg,egr Qg,exh ω eng ] .The controller design objective can now be formulated as the synthesis of a control algorithm that maximizes the WHR system output power while guaranteeing safe operation, i.e., χ f,mix ≥ 1, in the presence of highly dynamic engine disturbances.However, an optimization problem that maximizes power with χ f,mix ≥ 1 gives a non-unique global solution as a function of control input u.This could lead to excessive system actuation [21].Additionally, in order to avoid infeasibility during optimization, the output constraints need to be included as soft constraints, which requires additional optimization variables.Thus, the problem complexity increases.To avoid these undesired effects while satisfying the safe operation requirement, here, we use a particular solution, in which the vapor quality after both evaporators indicates vapor, i.e., [χ f,egr χ f,exh ] ≥ [ 1 1] .The expander is able to cope with working fluid in the two-phase state.Thus, for shorts periods of time, we assume a vapor fraction limit of 0.9 at the expander inlet.The WHR system efficiency increases as the system operates closer to the safety margin [10], i.e., χ f,egr = χ f,exh = χ f,mix = 1.Consequently, based on simulation results, the WHR output power is shown to reach its maximum at the safety margin [21].Thus, for maximum output power, the controller design objective reduces to the synthesis of a control algorithm with good disturbance rejection properties that allows operation close to the safety limit.

Model for Control Design
For control design, a low-computational model that captures the WHR system dynamics is necessary.However, the developed high-fidelity WHR system model has around 150 states and is highly nonlinear, mainly due to the two-phase phenomena.Therefore, the model is linearized and reduced around a representative set of operating conditions.In Figure 6a, we show the steady-state WHR system operating area, including limitations: the sum of Qg,egr and Qg,exh is limited to 100% by the condenser cooling power; ethanol decomposition can occur at a temperature higher than 280 • C, and emissions can increase due to the low exhaust gas flow through the EGR valve.These limitations are expressed in the steady state, and thus, they can be violated for short periods of time, especially during transients.The resulting WHR system operating area is shown in green in Figure 6a.
Within the WHR operating area, three linearization points are selected (see Figure 6a with numbers) that are shown to be sufficient in describing the WHR system dynamics [16].The points are selected close to the boundaries for emissions and condenser limitation, where the system is expected to operate more often.This is seen from a World Harmonized Transient Cycle and steady-state measurements that also motivate the switching lines choice l 1 and l 2 .
The three corresponding linear models are of high order, because they retain the same number of states as the high-fidelity model.To use these models for control, model reduction techniques based on balanced truncation are used.A model order of n = 7 is chosen, such that the error system ||G − G r || ∞ satisfies: where G is the full linear model of order m, G r is the reduced linear model of order n and σ are the singular values of G.The reduced order state-space linearized model in discrete time is: where The three reduced order models are integrated within the control scheme using a switching mechanism with hysteresis, as shown in Figure 6b.The switching mechanism indicates the active controller based on the EGR and EXH heat flow rates from the engine.This is achieved by monitoring the distance between the switching lines l 1 , l 2 and the point that expresses the current EGR and EXH heat flow rate.A parameter δ q that defines the hysteresis is introduced to avoid chattering around the switching lines.From simulation results, we conclude that δ q = 14.5% gives good switching behavior without chattering.

Linear Model Predictive Control
The control inputs u(k) = [u egr (k) u exh (k)] , at time k, are computed by an MPC controller.The controller consists of a state estimator and an optimizer (see Figure 5).The state estimator is used to obtain estimates of x(k) from the output y(k) = [χ f,egr (k) χ f,exh (k)] computed from the measurements.To improve the state estimation, the engine disturbances w(k) are assumed to be measured.
The optimizer solves a linear optimization problem that minimizes a cost function J subject to constraints.As the objective is to have good disturbance rejection, the cost function J penalizes deviations of the predicted controlled outputs ŷ(k + i) from a constant reference r(k + i) and changes of the inputs ∆u(k + i).The penalty matrices W y and W ∆u are constant inside each operating region of the WHR system.The predictions are expressed in terms of ∆u(k + i) = u(k + i) − u(k + i − 1), and the cost function is defined as: where N y is the prediction horizon and N u is the control horizon.The cost function J(k) is obtained based on the reduced order linear model, which is used to predict the behavior of the plant, starting at the current time k, over a future prediction horizon N y .This predicted behavior depends on the assumed input trajectory u(k + i) over a control horizon N u and measured output y(k) and disturbance w(k).
The input increments ∆u(k In Equation (30), || • || 2 W denotes the weighted two-norm squared, i.e., ||x|| 2 W = x W W x. At each time k ∈ Z + , we estimate the system state x(k) and solve the following optimization problem (this is often called a receding horizon principle): subject to Equation (29) and: The control action sent to the system is u(k) = u(k − 1) + ∆u * (k), where ∆u * (k) is the first element of the optimal sequence obtained by solving Problem 1.For state estimation, a general linear state observer [22] is designed based on the measured plant output, which gives the state estimate x through a gain matrix.The gain matrix is designed using well-known Kalman filtering techniques.
For each region in the operating WHR area, an MPC controller is designed.All three MPC controllers receive the current manipulated variable and output signals, for their state estimates.During switching, the state estimator can experience bumps, because the model reduction does not preserve the same states for each region.However, by penalizing the input increments ∆u(k +i|k) in the optimization Problem 1, these effects are mitigated and aggressive plat actuation is avoided.

Nonlinear Model Predictive Control
In this section, a nonlinear model predictive control (NMPC) is developed using the high-fidelity WHR system model.First, let us denote the high-fidelity WHR system using the general form of a non-linear system: where x is the state vector, u the control input, w the disturbance vector and y the model output.Then, at each time step k ∈ Z + , the NMPC strategy solves the following optimization problem: subject to Equation (33) and: The NMPC uses the same prediction horizon N y and control horizon N u as the linear MPC.However, note that Problem 2 computes directly the control input u(k) over the control horizon N u and not the control input increment ∆u(k).Moreover, a state estimator and a switching mechanism are not necessary, since the control model and the high-fidelity model are identical.
The reason for developing an NMPC strategy is only for comparison purposes, since the algorithm is computationally expensive and, thus, is not feasible for online implementation.

Simulation Results and Discussion
To evaluate the MPC strategy, measurements from a Euro VI heavy-duty diesel engine are used.The control strategy is tested on a simple stepwise cycle and then on a cold-start World Harmonized Transient Cycle (WHTC) that covers typical driving conditions: urban, rural and highway.As a benchmark, a classical PI controller is used [16], as well as a nonlinear MPC (NMPC) strategy.The PI control objective is to manipulate the control input u, such that the evaporator's outlet temperature follows a predefined reference temperature.This reference temperature is the ethanol saturation temperature plus a safety margin of 10 • C. The simulation environment is MATLAB/Simulink ® from The Mathworks Inc., Massachusetts, MA, USA.The MPC strategy is implemented using the model predictive control Toolbox, while the NMPC strategy uses a sequential quadratic programming (SQP) solver.
The controller sampling time is T s = 0.4 s, which yields an acceptable approximation of the system dynamics.The prediction horizon is N y = 45, and the control horizon is N u = 4 time steps.The prediction horizon of 18 s is chosen, so that the most important system dynamics within the engine operating area are captured.The control horizon of 1.6 s follows from the fact that larger values do not considerably improve the performance, while the computational complexity increases significantly.The corresponding weighting matrices are: MPC {1,2,3} : W y = {diag(0.8,0.8), diag(0.9,0.9), diag(0.8,0.9)}, W ∆u = {diag(0.05,0.05), diag(0.1,0.1), diag(0.2,0.2)}, NMPC: W y = diag(0.9,0.8), For the MPC strategy, the weighting matrices are tuned for a stepwise cycle, such that a good trade-off is obtained between control effort and robustness to model uncertainty.However, using these settings, the controller is shown to perform well also for other type of cycles, such as a WHTC.For the NMPC strategy, a similar weighting is chosen for W y .Since the NMPC optimization is performed subject to the high-fidelity model in Equation (33), the weighting W u is chosen small to allow more control freedom.

Stepwise Cycle
In this cycle, illustrated in Figure 7a, the engine torque is varied stepwise within the complete operating area, i.e., 500-2600 N•m, while the speed is kept constant at 1213 rpm.From the WHR system perspective, the engine torque and speed translates into EGR and EXH heat flow rate (see Figure 7b).The EGR and EXH heat flow rate together with the engine speed are the main disturbances for the WHR system, besides other disturbances, such as ambient temperature and ethanol temperature from the reservoir.The operating region is determined based on the EGR and EXH heat flow rate outside the controller.All three controllers become active with no chattering during switching, which could lead to instability.This indicates a good switching behavior of the switching mechanism with hysteresis.
In Figure 8, we show the vapor fraction after the EGR and EXH evaporator, as well as after the mixing junction.The reference r for the controller is set to a constant value of 1.05 superheat for both evaporators.From Equation (13), the superheat in temperature corresponding to this vapor fraction varies as a function of pressure from 10 • C to less than 5 • C as the pressure reaches its maximum.After the mixing junction, a good disturbance rejection is obtained, with ±5% maximum deviation from the specified reference r.Starting from 2400 s, the control inputs reached the limits (see the lower plot in Figure 8), i.e., u egr = 0% (fully closed) and u exh = 100% (fully opened).Thus, the output is indeed no longer able to follow the reference.The reason comes from the engine heat flow rate: at the EGR evaporator, the exhaust gas heat flow rate reached more than 60% while at the EXH evaporator, it dropped close to zero due to u g 2 actuation to avoid exceeding the condenser cooling capacity.However, depending on the situation, the vapor fraction after the mixing junction can still show the vapor state, and therefore, mechanical power is delivered, despite the control input that hit the constraint.

World Harmonized Transient Cycle
To validate and quantify the performance of the developed control strategy, a highly dynamic cold-start World Harmonized Transient Cycle (WHTC) is used (Figure 9a).This is an 1800 s cycle that consists of urban, rural and highway driving conditions.The WHTC is a challenging cycle due to a highly dynamic behavior that includes also the heat-up phase of the engine and WHR system.The engine disturbances from such a cycle are illustrated in Figure 9b.The switching signal indicates that only Regions 1 and 2 are reached.For comparison, we use two additional control strategies: a classical PI control scheme [16] and an NMPC strategy described in Section 4. The results for these three control strategies are shown in Figure 10.For brevity, we only show the vapor fraction after the mixing junction.
During the first 1200 s, the WHR system is heating up, and therefore, no vapor state is encountered after the mixing junction.Note, however, that for each evaporator, the heat-up period is different, up to 600 s for the EGR evaporator and 1200 s for the EXH evaporator.The exhaust evaporator takes longer to heat-up due to more mass and the aftertreatment system located upstream.The difference in heat-up period can be observed also from the period of time the control input u egr and u exh are on the constraint.
The objective of all three controllers is to maintain vapor and to be as close as possible to χ f,egr = χ f,exh = 1 for maximum output power.However, due to highly dynamic disturbances and the limitations of the control input, maintaining the vapor state is challenging.The main reasons for introducing an NMPC controller are: first, to verify how good is the control strategy if the WHR system model is exactly known, and second, what the global control input solution is.From the simulation results, we found that the global solution is difficult to obtain, since the system is highly nonlinear with local minima over the prediction horizon.Thus, the NMPC solution is sensitive to the initial condition for the optimization Problem 2. Here, for initialization of the NMPC strategy, we used the MPC solution.As a result, the NMPC strategy corrects for the modeling errors between the reduced order linear model and the high-fidelity model.At first instance, all three controllers seem to behave quite similarly, especially the MPC and NMPC control strategies.To quantify the performance of each controller, we use two indicators: the time in the vapor state t v and the recovered thermal energy U over the complete cycle.
The performance indicators are define as follows: Time in vapor state: Thermal energy: U = where t f is the total cycle time and P whr (t) is the instantaneous WHR system net output power.The instantaneous output power P whr (t) and the recovered energy U are calculated only when the safety requirement χ mix ≥ 1 is satisfied.In the case that the safety requirement is not fulfilled, i.e., χ mix < 1, the expander needs to be bypassed to avoid damage.
From Figure 11, we can see that the MPC controller outperforms the PI controller in terms of time in the vapor state and recovered thermal energy, with approximately 15%.The NMPC strategy improves the considered performance indices by 10% as compared to the MPC strategy.This result gives an indication about the model uncertainty used within the MPC strategy.Despite the improved result, the NMPC controller has several drawbacks.First, the computational complexity is significant, which makes the NMPC scheme not feasible for real-time implementation.Second, it can converge to local minima, resulting in worse control behavior than with a linear MPC.As a solution, an NMPC algorithm that uses global solvers can be considered; however, the computational complexity will grow even more.

Conclusions and Future Research
The modeling and control of a waste heat recovery system for a Euro-VI heavy-duty truck engine was considered.The exhaust gas energy is recovered by means of two parallel evaporators.To generate power, an expander is used, mechanically coupled with the pumps to the engine crankshaft.The main conclusions of this paper are summarized in what follows.
(1) The existing WHR system model was improved by including enhanced evaporator models.The model combines the finite difference modeling approach with a moving boundary one and, thus, captures multiple phase transitions along a single pipe flow.(2) A switching model predictive control strategy was designed, to guarantee the safe operation of the WHR system.The proposed control strategy performance is demonstrated on a highly dynamic cold-start World Harmonized Transient Cycle.Simulation results showed improved performance for the proposed MPC strategy in terms of vapor time and recovered thermal energy by 15%, as compared to a classical PI controller.(3) A nonlinear model predictive control was developed to demonstrate that control performance can be improved by approximately 10% in the case that the model of the system is exactly known.(4) One limitation of the proposed MPC strategy is the need of vapor fraction measurement equipment or an estimator to make the method applicable in practice.
Future work will focus on vapor fraction estimation and on experimental demonstration of the proposed MPC strategy.Furthermore, the supervisory control strategy [23] that combines energy and emission management will be extended by considering the proposed WHR control strategy.

Figure 1 .
Figure 1.Schematic representation of the experimental set-up.

Figure 2 .
Figure 2. Ethanol temperature as a function of specific enthalpy and pressure.

Figure 3 .
Figure 3. (a) Total exhaust gas heat flow rate, including losses and working fluid heat flow rate for both evaporators; (b) Working fluid mass flow rate through the evaporators and system pressure.

Figure 4 .
Figure 4. (a) Working fluid temperature after the evaporators and mixing junction; (b) Vapor fraction after the evaporators and mixing junction.

Figure 5 .
Figure 5. Schematic representation of the control scheme.

Figure 6 .
Figure 6.(a) Cold-start World Harmonized Transient Cycle and steady-state points within the waste heat recovery (WHR) system operating area (green) and limitations: cooling capacity (grey), high emissions (yellow), ethanol decomposition (magenta); (b) Schematic representation of the switching mechanism with hysteresis.

Figure 7 .Figure 8 .
Figure 7. (a) Engine torque and speed; (b) Engine exhaust gas heat flow rate and the operating region signal.

Figure 9 .
Figure 9. (a) Engine torque and speed from a cold-start World Harmonized Transient Cycle (WHTC); (b) Engine heat flow rate from a cold-start WHTC and operating region signal.

Figure 10 .
Figure 10.(a) Vapor fraction after the mixing junction, absolute error and control inputs for proportional-integral (PI), model predictive control (MPC) and nonlinear MPC (NMPC) control strategy for a cold-start WHTC; (b) WHR system net output power, system pressure and mass flow rates through the EGR and EXH evaporator.

Figure 11 .
Figure 11.Time in the vapor state and recovered thermal energy for the MPC and NMPC strategy compared to the PI control strategy.