mCHP Optimization by Dynamic Programming and Mixed Integer Linear by chenmeixiu


									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
 Stirling engine commercialized for example by Whispergen,                                                           3-way
                                                                                   Fuel valve
 Microgen and Solo. Stirling engines are more reliable and
 less noisy than internal combustion engines. Fuel cells will                                                 Ti
 be the next products to appear on the mCHP market with                            Fuel valve
                                                                                                                                                Hot Water
 even less maintenance costs and more promising electrical                                                                                       Demand
 efficiency (35 to 50 %).                                                                       Electricity
    With the liberalization of the Gaz and Electricity market,                                                                               Mixing
 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@                                          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@
    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:                                                although inherent problems are obvious. The demands are

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.
                                                                                                                 Open Loop
                            Weather Forecast                                                                                                    Open loop
                                                                                                                                     Model      trajectories
                            Needs                                                              T                       Open loop
                          Prediction                                                                                    Control
         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
                                                                                 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:
 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
                                                                                       *                                     *
 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
                                                                                      *                                  *

 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

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
                                                                                 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

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 )]
                                                                                                  ∆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
          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.

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
                                                                                                                      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                                   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
                                                                 Dynamic Programming                                                                   Dynamic Programming
                                                 80                                                                                       80
                     Input Temperature Ti (°C)

                                                                                                               Tank Temperature Tt (°C)

                                                 20                                                                                       60
                                                      0      4       8      12      16         20   24                                         0   4      8      12      16       20      24

                                                                  Linear Programming                                                                    Linear Programming
                                                 80                                                                                       80
                                                 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)

                                                                                                               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
                                                 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.

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)
                                      20                                                                            60
                                           0   4       8      12   16    20   24                                         0    4        8      12       16      20       24

                                                    Linear Programming                                                             Linear Programming
                                      80                                                                            80
                                      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)


                                                                                         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
                                      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.


To top