VIEWS: 52 PAGES: 6 POSTED ON: 6/30/2011 Public Domain
The 14th International Conference on Intelligent System Applications to Power Systems, ISAP 2007 November 4 - 8, 2007, Kaohsiung, Taiwan mCHP Optimization by Dynamic Programming and Mixed Integer Linear Programming D. Faille, C. Mondon, and B. Al-Nasrawi Abstract—This paper studies the optimal control of a optimization of a district was considered. The system domestic micro combined heat and power system. The objective consisted of 200 houses and buildings, a fleet of mCHP and function is to reduce the gas and electricity bill of the storage devices. The commitment of the machines and the consumer. The problem is nonlinear and contains continuous storage strategy was solved by Mixed Integer Linear and logical variables. Two solutions are compared in this Programming (MILP). In [4], the problem of the local paper: Dynamic Programming and Mixed Integer Linear Programming. optimization of a mCHP was addressed. A predictive control was proposed to achieve a compromise between the cost I. INTRODUCTION minimization and the thermal fatigue due to start and shutdown transients. The solution was given by Dynamic Micro Combined Heat and Power Plant (mCHP) are new Programming (DP). In the present paper, we go further with devices that produce Heat and Electricity for single and this latter problem by proposing a solution with MILP multifamily houses or office buildings and that could in the solver. In the second section, we present the system we want future replace the traditional boilers to fulfill the heating and to control, the general structure of the proposed control the hot water requirements as well as a part of the electricity scheme and a brief overview of the two algorithms DP and demand (the connection to the grid for a dwelling equipped MILP. The third section presents the modeling. Simulation with mCHP in most cases remains necessary). results and comparisons are given in the last section before a The mCHP solution that is being developed today in conclusion that indicates the direction for future work. Europe (e.g. United Kingdom, Germany and France) presents several benefits. First, it contributes to the II. GENERAL SETUP Greenhouse Gas Reduction thanks to a better global efficiency. Second, it diminishes the consumer’s energy bill A. The Process with savings that can reach in certain cases 30% [1],[2]. The system that we consider is shown in Fig. 1. It consists During certain periods when the local demand is low, the of a house equipped with a mCHP, an auxiliary boiler, a electricity produced by the mCHP can be sold to the central heating system and a hot water cylinder. The mCHP Distribution Network Operator and thus give income to the and the auxiliary boiler are used to cover the heating and the consumer. hot water requirements. The mCHP produces electricity that Different mCHP technologies are present on the market can be consumed locally or exported. For more details about today. The first is the gas engine technology with vendors the process, the reader can refer to [4]. such as Honda, Ecopower, Senertec, etc. This technology is well proved but the maintenance costs are still high. A more Thermostatic To Valve recent technology is the external combustion engine or Auxiliary Heat Demand Boiler Stirling engine commercialized for example by Whispergen, 3-way valve Fuel valve Microgen and Solo. Stirling engines are more reliable and mCHP less noisy than internal combustion engines. Fuel cells will Ti be the next products to appear on the mCHP market with Fuel valve Tt Hot Water even less maintenance costs and more promising electrical Demand Cylinder efficiency (35 to 50 %). Electricity Demand With the liberalization of the Gaz and Electricity market, Mixing valve the French Electricity Company EdF is interested in solutions that combine these two energies. EdF R&D has Fig. 1. Plant Description launched several projects on these topics. In [3], the B. The control The general scheme we propose in order to optimize the D. Faille is with the Electricité de France Research and Development, 6 Quai Watier, 78400 Chatou, France (phone: 33-1-30-87-84-63; fax 33-1- process is given in Fig. 2. It consists of two modules. 4087-; e-mail: damien.faille@ edf.fr). The first module calculates the forecasts for the heating, the C. Mondon, with the Electricité de France Research and Development, 6 hot water and the electricity demands. These demands Quai Watier, 78400 Chatou, France (e-mail: damien.faille@ edf.fr). Boris Al-Nasrawi is with the European Institute for Energy Research depend mainly on the weather and the consumers behavior. (EIFER), Emmy Noether Strasse 11, 76131 Karlsruhe, Deutschland , (e- This module is beyond the scope of the present paper mail: boris.al-nasrawi@ edf.fr). although inherent problems are obvious. The demands are 572 The 14th International Conference on Intelligent System Applications to Power Systems, ISAP 2007 November 4 - 8, 2007, Kaohsiung, Taiwan difficult to model particularly for the electricity demands Prices, Needs that are characterized by peaks. Xo(k) Open Loop Weather Forecast Open loop Optimization Model trajectories Needs T Open loop Prediction Control Feedback Needs Forecast Gas and Electricity Prices Control Optimal State Perturbation Estimator Process Control Perturbations - Weather - Inhabitant Fig 3. Predictive Model based Control. mCHP Heat +Support Boiler + Hot water C. Optimization Background +Control +Electricity 1) Dynamic Programming (DP). DP is a classic +Control mathematical method which aims to find step by step the PROCESS best decisions to control a dynamic system. Indeed, it is a Fig 2. General Control Structure way to solve the tradeoff between immediate low present The second module that is treated in the present paper is cost and undesirable high future cost. The framework of DP the optimal control of the process. The control policy was defined by Bellman [7], and recently Bertsekas wrote a calculates the control inputs i.e. the opening of the 3-way complete state-of-the-art book on this technique [8]. The valve, the loads of the mCHP and of the auxiliary boiler, so horizon can be finite or infinite; the problem can be that the needs are satisfied and that the bill is minimized. deterministic or stochastic; time, state and control can be The optimization horizon is finite and in this study is 24 discrete or continuous. However, in every case, the system hours. under consideration must have two mandatory properties: The optimal control must fulfill some practical additive over time cost function, and a first order dynamic requirements. The first requirement is related to the transition equation. In our study, we consider a finite computation time. With the development of domotics discrete time horizon and a deterministic situation. Hence system, the homeowners have at their disposal increasingly the problem can be formulated as follows: N powerful computers to control home applications, however min ∑ C (t , xt , ut ) + C F ( x N ) the solution of the optimization problem must be found πk t =k within reasonable delay. This requirement limits therefore xt +1 = f (t , xt , ut ) , ∀ t ∈ { k , … , N − 1 } the complexity of the optimization algorithm that can be Where t is the discrete time index, xt and ut are used. respectively the state and control variables, f is the transition The second feature that the control must respect is the law and πk ={µk,…, µt,….µN } is a control policy. Indeed, we robustness. The optimization is based on several models. are not only looking for optimal trajectories but also for the The process model describes the dynamic behavior of the * controlled outputs in response to given control inputs and is optimal state feedback function uk =µk(xk), which relates the always an approximation of the reality. The perturbation control input to the state of every sample k. models that correspond to the Needs Prediction module in If π k ( x k ) is the optimal solution for this problem and * Fig. 2 are difficult to obtain and are prone to errors too. * * Two techniques can be adopted to mitigate the impact of the ν k ( x k ) = J (π k ( x k )) the Bellman Value, the principle of errors on the optimization. The first one is the receding optimality leads to the following recursive relation called the horizon principle (see [5] for a presentation) that is a key Bellman equation: stone of the predictive control theory. The idea is to find a ν k ( x k ) = min C (k , x k , u k ) + ν k +1 ( f (k , x k , u k )), ∀x k * * uk sequence of control inputs that optimizes the criterion over a finite horizon and apply only the first point. Then we and estimate the state and launch a new optimization with a u k = arg min C (k , x k , u k ) + ν k +1 ( f (k , x k , u k )), ∀x k * * uk shifted time window. The state estimator is represented in the Fig. 3. The second technique is proposed in [6] and is Starting from the end (k =N), it is possible to compute the integrated in our control scheme shown in Fig. 3 too. The Bellman Values and the optimal control step-by-step for control inputs calculated by the optimization are applied to a every state and finally find the solution for the first-step. model to give the reference trajectories that are compared to Practically, in most cases, state and control variables are the real trajectories. The difference is reduced by a feedback discretized, thus Bellman equation can be solved by simple controller (PID, relay, etc.) that calculates a correction term enumeration. added to the optimization results. 2) Linear Programming (LP) and Mixed Integer Linear Programming (MILP). LP is another popular optimization 573 The 14th International Conference on Intelligent System Applications to Power Systems, ISAP 2007 November 4 - 8, 2007, Kaohsiung, Taiwan method, which has been introduced by Dantzig and his continuous solutions. Simplex algorithm [9]. The framework of LP deals with any Considering now computation time, DP suffers from the so- optimization problem, under the condition that every part is called curse of dimensionality. Indeed, complexity of DP is composed of linear terms: in O(Card(S).Card(C).Card(T)) where S, C and T are min c t .x respectively the state, control and time set and Card(X) is the x cardinal of the set X. Therefore the computation time grows A.x = a exponentially with the number of control and state variables. B.x ≤ b In practice, three variables is an acceptable size. Research Where x is the vector of continuous optimization has been done for many years to overcome this difficulty, unknowns, c is the cost vector, A and B are matrices of with techniques like Lagrangian relaxation [11] or equalities and inequalities. Indeed, A and B define a convex Approximate Dynamic Programming (see [8], pp 283-386). polytope in the variable space and –c is a constant Incredibly huge LP problems can be solved and has been minimization direction. Then it can be proved that the proved to have polynomial worst-case complexity with new optimal solution lies in one of the vertices of the polytope. algorithms: ellipsoid method [12], or interior point method The Simplex algorithm explores the set of vertices with [13]. MILP is NP-hard: the binary tree has an exponential simple and fast computations, by a greedy strategy to jump size, and in most practical case a good solution rather than from one vertex to an adjacent one with a local maximum the global optimum is computed. By using cutting plane minimization. Even, if in theory the whole set might be techniques, constraint propagation and even heuristics, explored, in practice the Simplex algorithm finds the today’s commercial solvers have made impressive progress optimum without exploring all the vertices. and are now able to solve in a reasonable time many MILP is an extension of LP where a part of the decision problems [10]. variables are binary. In this way, it is possible to take into To conclude this comparison, whereas MILP find one account discrete values or to express logical conditions. solution or trajectory in optimal control problems, DP Basically, the methods applied to solve MILP problems in computes a control law for each time step. Hence, DP can practice use a two-level algorithm to find a solution. At the correct modeling errors or perturbations with no extra- upper level, the tree of the binary variables is explored by a computation, so information about future is unchanged. Branch & Bound method. At each step of this method, a part of the binary variables are set, and a continuous problem is III. MODELING formulated by relaxing the integrity condition of the other In this section, we propose a modeling of our control free binary variables. A Simplex routine is used to solve this problem in the DP and MILP framework. Basically the relaxed problem. Working with the lower and upper bound, models contain a cost function that is corresponding to the given respectively by the relaxed problem and the best electricity and gas bill and a state equation that corresponds solution currently found, only a part of the tree is explored to dynamic behavior of the process. (see [10]). 3) DP vs. MILP. Both DP and MILP can be used to find A. Modeling for Dynamic Programming. an optimal control of the mCHP. Here is a brief comparison The system described in Fig. 1 can be modeled by detailed of the two methods to show their advantages and their models based on physical equations. However for control drawbacks. purpose, the model must achieve a compromise between First of all, the problem must fit the framework of each accuracy and tractability, and its complexity must be kept as method. The problem modeling can be more or less complex low as possible. according to the method used. For instance because of the The mCHP is supposed to be a combustion engine that can’t first order equation assumption the DP problem formulation adjust its load. Its thermal and electrical efficiencies sometimes necessitates the use of state augmentation η thm , η em can therefore be constant. They are defined as the techniques (see [8] p 35), which increase the computation ratio between the electrical and the thermal power produced time. With MILP, the linearity assumption is often and the natural gas power consumed by the machine Qgm . restrictive. One dimension nonlinearity can be well approximated by piecewise linear functions but The auxiliary boiler is supposed to adjust its load and is nonseparable nonlinearity such as product are much more characterized by a thermal efficiency constant η thsb . Its gas difficult to handle. In every case, approximations used to power consumption is noted Qgsb . overcome the nonlinearity result in an increased number of The thermal power produced by the mCHP and/or the binary variables and a consequent increase in computation auxiliary boiler are either used to fulfill the heating demand time. Conversely, DP handles naturally every nonlinearity Pd h or the hot water demand Pd hw . The repartition is done and MILP has no problem with time-coupling constraints. Furthermore, when solving the problem, a key issue in DP by the 3-way valve opening α k that is binary. The heating may be the discretization of continuous variables, which lead and hot water temperature responses are not instantaneous, to an approximation, whereas MILP computes exact but are influenced by the inertia of the circuits: Ct is the 574 The 14th International Conference on Intelligent System Applications to Power Systems, ISAP 2007 November 4 - 8, 2007, Kaohsiung, Taiwan thermal inertia coefficient of the hot water circuit and Pesk + Pebk = Pdek + n em Q gmk (9) depends mainly on the cylinder design, Ci is the thermal Now let define Pth and Pthw the part of the total thermal inertia coefficient of the heating circuit. The inertia power used for the heating and for the hot water. coefficients are defined as the ratio between the energy Pt h k = (1 − α k )(η thm Q gmk + η thsb .Q gsbk ) (10) increase during a time step ∆t and the resulting temperature Pt hwk = α k ⋅ (η thm Q gmk + η thsb .Q gsbk ) (11) increase of the water at the input of the mCHP for the heating circuit ( ∆Ti ) and in the hot water tank for the hot These new variables are defined as a product of a logical water circuit ( ∆Tt ). variable αk and a continuous term ηthm Qgmk + ηthsb .Qgsbk that The system described in Fig. 1 can thus be modeled by the is the total thermal power produced. This quantity takes its following set of equations. Equations (1-2) represent the values between a maximum M and a minimum m. With this dynamic responses of the heating and the hot water circuits. assumption, nonlinear expressions (10) and (11) can be Equations (3-4) are the corresponding temperature limits. replaced by linear inequalities (12-15) and (16-19). Equations (5-7) are the constraints on the control inputs. Ttk +1 − Ttk 1 = [− Pdhwk + Pthwk ] (1b) Tt k +1 − Tt k 1 ∆t Ct = [− Pd hwk + (1) ∆t Ct Tik +1 − Tik 1 = [− Pdhk + Pthk ] α k ⋅ (η thm Q gmk + η thsb .Q gsbk )] (2b) ∆t Ci Tik +1 − Tik 1 Pthhk ≤ (1 − α k )M (12) = [− Pd hk + (2) ∆t Ci Pthhk ≥ (1 − α k )m (13) (1 − α k ) ⋅ (ηthm Qgmk + ηthsb .Qgsbk ) ] Pthhk ≤ ηthm Qgmk + ηthsb .Qgsbk − mα k (14) Ti min ≤ Ti k ≤ Ti max (3) Pthhk ≥ ηthm Qgmk + ηthsb .Qgsbk − Mα k (15) Tt min ≤ Tt k ≤ Tt max (4) Pthhwk ≤ α k M (16) α k ∈ { 0 ,1} (5) Pthhwk ≥ α k m (17) { Qgmk ∈ 0 , Qgm max } (6) Pth hwk ≤ η thm Q gmk + η thsb .Q gsbk − m(1 − α k ) (18) 0 ≤ Qgsbk ≤ Qgsb max (7) Pth hwk ≥ η thm Q gmk + η thsb .Q gsbk − M (1 − α k ) (19) The sample time is in hour. The temperatures are in °C, the gas power in kW, the inertia in kWh/°C. C. Simplification It remains to define the objective function J in (8) that is If we suppose that the 3-way valve opening can be the bill of the consumer. We suppose that we know the continuously adjusted like the auxiliary boiler load can be, prices of the gas cg and of the electricity ce. We define a the problem is simpler for the DP as well as for the MILP. repurchase ratio τ as the ratio of the price of the electricity This simplification is valid only with large time step bought by the utility to the price of the electricity sold by the (hourly), otherwise a modeling error is introduced. utility. The prices are given in €/kWh. 1) MILP simplification. We keep (1b, 2b, 3, 4, 6, 7, 8b, 9), and we replace the set of inequalities (12-19) by the [c N −1 J = ∆t ∑ (Q gmk + Q gsbk ) + (8) k =0 gk following constraints on thermal power (20-22). { cek max( Pdek − η emQgmk ,0) + τ min( Pdek − ηemQgmk ,0) } ] Pthk + Pthwk = η thm Q gmk + η thsb .Q gsbk (20) The goal of the control is to find the load of the machines Pthk ≥ 0 (21) Qgmk, Qsgbk and the opening of the valve αk that minimize the Pthwk ≥ 0 (22) bill (8) and respect the constraints (1-7). This problem can The variable αk has been eliminated but we can be solved by DP once the continuous variables (Ti, Tt and reconstruct it knowing Pthk and Pthwk when the total power Qgsb ) have been spatially discretized. is non zero. Thanks to this simplification, a binary variable B. Model Formulation for MILP solution. is suppressed and the computation time much faster. To be able to solve the problem by MILP we have to get 2) DP Simplification. When we consider that αk and Qgsbk rid of the nonlinearity. The first one is the product between can be continuously adjusted, we can find an exact solution the inputs in (1) and (2). The second is the min and max for the quantities Pthk and Pthwk that respect (1b) and (2b) for term in (8). This is done by introducing auxiliary variables. every Tt and Ti defined on the grid of discretization at time k Let define Pes and Peb the sold and bought electricity and k+1. So the DP algorithm doesn’t require any power, (8) becomes (8b). interpolation, because the trajectories are going through points of discretization. Numerical errors are reduced, and a [ N −1 J = ∆t ∑ c gk (Q gmk + Q gsbk ) + c ek ( Pebk + τ Pesk ) ] (8b) coarser discretization grid can be used, implying faster k =0 The new variables are defined by the following relation. computation time. 575 The 14th International Conference on Intelligent System Applications to Power Systems, ISAP 2007 November 4 - 8, 2007, Kaohsiung, Taiwan IV. RESULTS constrained between 60°C and 80°C. The mCHP and In this section, we consider the following application. The auxiliary boiler efficiencies are ηthm= 0.65, ηem=0.2, and demands correspond to a mid-season day. We can see their ηthsb=0.82. The auxiliary boiler gas power consumption evolutions on the curves given in Fig. 4. The time sample is varies between 0 and 20 kW. The mCHP gas power 10 minutes. consumption has two levels 0 and 5 kW. For DP, the temperatures have been discretized with a pace of 1°C, and the auxiliary boiler consumption has been discretized with a electricity pace of 0.1 kW. Demand (kW) 10 heating The DP algorithm has been coded in Matlab™, and the hot water MILP problem has been coded in the software XPRESS- MP™. The two problems have been solved on the same 5 computer. The computation time is 6.6 s for DP and 63.3 s for MILP. The trajectories obtained with the two methods given in 0 0 4 8 12 16 20 24 Fig. 5, are very similar. The cost function is 3€26 for the DP Time (h) solution and 3€25 for the MILP solution. Fig 4. Electrical , heating and hot water demands. We observe in the Fig. 5 that the auxiliary boiler is not used The price of electricity ce has two levels. The high price is on the transient, thanks to the energy storage and destorage 11 c€/kWh between 7 am and 10 pm; the low price is 5 in the heating circuit and in the tank. We observe that the c€/kWh. The repurchase ratio τ is 4:11. The gas price cg mCHP transients contain several startups and shutdowns. A equals to 2 c€/kWh. The heating circuit inertia Ci is equal to solution proposed in [4] to reduce the number of these 0.06 kW/°C and the hot water circuit inertia Ct is equal to transients is to add a penalty term to the criterion. Another 0.46 kWh/°C. These inertia coefficients have been estimated possible approach is to add a constraint on the minimal with step responses on a detailed model [4]. duration of the on-position and off-position of the mCHP. A The temperature at the mCHP entrance is limited between simpler solution adopted in the present paper is to take a 20°C and 80°C. Similarly, the tank temperature is larger sample time, so the control inputs are changing less often. Dynamic Programming Dynamic Programming 80 80 Input Temperature Ti (°C) Tank Temperature Tt (°C) 60 70 40 20 60 0 4 8 12 16 20 24 0 4 8 12 16 20 24 Linear Programming Linear Programming 80 80 60 70 40 20 60 0 4 8 12 16 20 24 0 4 8 12 16 20 24 Time (h) Time (h) Dynamic Programming Dynamic Programming 6 20 mCHP gas power (kW) 4 boiler gas power (kW) 2 10 0 0 0 4 8 12 16 20 24 0 4 8 12 16 20 24 Linear Programming Linear Programming 6 20 4 2 10 0 0 0 4 8 12 16 20 24 0 4 8 12 16 20 24 Time (h) Time (h) Fig 5. Comparison of DP and MILP solutions with a 10-minute sample. 576 The 14th International Conference on Intelligent System Applications to Power Systems, ISAP 2007 November 4 - 8, 2007, Kaohsiung, Taiwan Dynamic Programming Dynamic Programming 80 80 Input Temperature Ti (°C) Tank Temperature Tt (°C) 60 70 40 20 60 0 4 8 12 16 20 24 0 4 8 12 16 20 24 Linear Programming Linear Programming 80 80 60 70 40 20 60 0 4 8 12 16 20 24 0 4 8 12 16 20 24 Time (h) Time (h) Dynamic Programming Dynamic Programming 6 20 mCHP gas power (kW) 4 boiler gas power (kW) 2 10 0 0 0 4 8 12 16 20 24 0 4 8 12 16 20 24 Linear Programming Linear Programming 6 20 4 2 10 0 0 0 4 8 12 16 20 24 0 4 8 12 16 20 24 Time (h) Time (h) Fig 6. Comparison of DP and MILP solutions with a 1-hour sample. The Fig. 6 shows the results of an optimization done with a larger sample time. In this simulation, we suppose that the REFERENCES auxiliary boiler and the opening of the 3-way valve can be [1] Entchev E. “ Residential fuel Cell energy systems performance adjusted between the minimum and maximum limits. optimization using soft computing techniques”, Journal of power Sources,vol. 118, pp. 212-217, 2003 Between two time steps, the feedback control presented in [2] Peacock A.D., Newborough M. “Impact of micro-CHP systems on Fig. 3 can compensate the perturbations. domestic sector CO2 emissions,” Applied Thermal Engineering, vol. The optimal cost obtained is 3€28 (DP) and 3€26 (MILP). 25, pp. 2653-2676, 2005 The transients are almost the same. The time calculation is [3] C . Mondon., D. Faille “Optimization of a micro Combined Heat and Power Fleet”, in Proc. IASTED, 2005. shorter for the MILP (0.1s) than for the DP (6s). We can [4] D. Faille, C. Mondon, L. Heinckes,”Optimization of a micro observe that the MILP solver is more sensitive to the time combined heat and power Plant”, in Proc. IFAC Power Plant and step size than the DP. Power System Conference, 2006 [5] A. Bemporad, M. Morari: “Control of Systems integrating logic, dynamics and constraints”, Automatica, vol. 35, pp. 407-427, 1999 V. CONCLUSION [6] K. Pitscheider, B. Meerbeck, E. Welfonder: “Robust Model-based unit control concept with regulated deactivation of preheaters and heat The paper presents optimizations of a mCHP done by condensers”, in Proc. IFAC Power Plant and Power System MILP and DP. The results are very similar. The DP is more Conference, 2000 suited to integrate nonlinearity and MILP can easily [7] R.E. Bellman, “Dynamic Programming”, Princeton, N J: Princeton University Press, 1957 integrate time-dependent constraints. This study shows that [8] D.Bertsekas, “Dynamic programming and optimal Control”, Athena DP compares well with MILP for our application. DP is scientific, 2005 simpler to be implemented in real-time application because [9] Dantzig, “Linear Programming and extension “, Princeton University Press, New Jersey, 1963. it doesn’t need a solver like XPRESS-MP™ and presents a [10] R. E. Bixby, M. Fenelon, Z. Gu, E. Rothberg, R. Wunderling., “MIP, constant computation time. Furthermore, stochastic features Theory and practice – closing the gap”, ILOG technical paper of thermal and electrical needs should be better handled in [11] A. Merlin, P. Sandrin, “A new method for unit commitment at Electricité de France”, IEEE Transactions on Power Apparatus DP framework. For these reasons, DP should be the selected Systems, 102(5), 1983, 1218-1225 technique for mCHP optimal control. [12] L. G. Khachian, “A polynomial algorithm in linear programming”, The future work is to implement and validate the control Doklady Akad. Nauk SSSR, vol. 244, No. 5, pp.1093-1096,1979. scheme on a real-time test facility at EIFER and develop a [13] I. Adler, N. Karmarkar, M. Resende, G.Veiga, ”An implementation of Karmarkar's algorithm for linear programming”,. Mathematical needs prediction module. Programming, vol. 44, pp. 297-335, 1989. 577