Titre: Conception optimale de systèmes multi-énergies locaux : modèles de programmation linéaire à nombre entier mixte et approches de décomposition à deux niveaux Mots clés: Système multi-énergie, Programmation linéaire en nombres entiers, Optimisation de planification Résumé: Un système multi-énergie local (LMES) est un système énergétique décentralisé produisant de l'énergie sous de multiples formes pour satisfaire les besoins énergétiques de clients situés dans son voisinage. Les LMES présentent de nombreux avantages, par exemple une efficacité de production plus élevée et des coûts de maintenance plus faibles. Toutefois, afin de fournir leurs meilleures performances potentielles, les LMES doivent être conçus avec soin.

Le problème de conception d'un LMES, en tant que problème d'optimisation combinatoire, est difficile à résoudre en raison de sa grande taille. De nombreux outils numériques pour aider les gens à concevoir des LMES reposent sur des hypothèses fortes et des simplifications, et ils limitent généralement la taille des instances considérées afin de pouvoir proposer des solutions dans un temps de calcul acceptable.

Cette thèse de doctorat se concentre sur la modélisation et la résolution des problèmes de conception optimale des LMES. En termes de modélisation du problème, nous permettons de construire le système en plusieurs phases d'investissement et de considérer l'opération détaillée aux pas de temps horaires. Le problème d'optimisation qui en résulte, sous la forme d'un grand programme linéaire à nombres entiers mixtes, est résolu par deux approches de décomposition : un algorithme de décomposition hiérarchique étendu et un algorithme de décomposition de Benders généralisé.

Les résultats numériques obtenus à partir de trois études de cas réels montrent d'abord que les algorithmes de décomposition proposés sont plus performants que les approches de résolution traditionnelles. Nos résultats montrent également que, même si certaines approximations sont appliquées, les plans de déploiement obtenus sont de très bonne qualité.

Title: Optimal design of local multi-energy systems: mixed-integer linear programming models and bi-level decomposition approaches Keywords: Multi-energy system, Mixed-integer linear programming, Planning optimisation Abstract: A local multi-energy system (LMES) is a decentralized energy system producing multiple forms of energy to satisfy the needs of customers in its vicinity. LMESs display many advantages, e.g., higher production efficiency and lower maintenance costs. However, to provide their best potential performance, LMESs should be carefully designed.

The design problem of an LMES, as a combinatorial optimization problem, is difficult to solve due to its large size. Many off-the-shelf numerical decision-aid tools rely on strong assumptions and simplifications, and they usually limit the size of the considered instances to enable the resolution process to terminate within an acceptable computation time.

This PhD thesis focuses on the modelling and resolution of optimal design problems of LMESs. In terms of problem modelling, we allow to construct the system within multiple investment phases and consider the detailed operation with hourly timestep. The resulting optimization problem, as a large mixed-integer linear program, are solved by two decomposition approaches: an extended hierarchical decomposition algorithm and a generalized Benders' decomposition algorithm.

The numerical results obtained from three reallife case studies first show that the proposed decomposition algorithms significantly outperform the traditional solving approaches. Our results also show that, even if some approximations are applied, the obtained deployment plans are of very good quality.

Acknowledgement

During the past three years of my doctoral study, I was fortunate to receive help and support from many people. This thesis could not have been realized without their help. I would like to express my sincere gratitude to them.

First, I would like to express my deepest gratitude to my supervisors of thesis: Dominique Quadri, Céline Gicquel, Côme Bissuel and François Courtot. Despite the difficulties caused by COVID-19, they helped me with patience and kindness, which allows me to carry out my work remotely. Especially the weekly meetings with Céline and Côme, the discussions and inspirations during these meetings are the key for me to accomplish this thesis.

I would also like to thank François Clautiaux and Pierre Fouilhoux for having accepted to be the reporters of my thesis. My gratitude also goes to Frédéric Babonneau, Safia Kedad-Sidhoum, Sandra Ngueveu and Guillaume Sandou for being part of the jury for my doctoral thesis. Special thanks to Wenhui Du, the director of EDF R&D China, for proposing this excellent opportunity to me and for supporting my work during my doctoral study. I would also like to thank my colleagues in EDF R&D and I appreciate the great time we had working together.

Finally, I would like to dedicate this thesis to my wife Xinran Wang, without whom I could not have accomplished this long journey. I would also like to thank my mother Linmei Zhang. Words cannot express my gratitude to her consistent encouragement and selfless love. I am also extremely grateful to my father Yuming Liu, who is the first and the best teacher in my life. I would also like to say to my two younger sisters, Qianqian Liu and Youyou Liu, I love you.

Abstract

Over the last decades, the energy industry has been striving to improve the energy production efficiency, to lower the greenhouse gas emissions related to energy production and distribution and to better integrate renewable energy resources. Local multi-energy systems (LMESs) are an interesting alternative to meet these challenging objectives. Basically, an LMES is a decentralized energy system producing energy within multiple forms to satisfy the energy needs of customers located in its vicinity. The customers correspond to a set of buildings belonging e.g. to a university campus, a hospital complex or a city district. LMESs have a higher production efficiency and a lower maintenance cost than traditional individual (i.e., singlebuilding) energy systems.

LMESs thus display many practical advantages. However, in order to provide their best potential performance, they should be carefully designed. Designing an LMES essentially consists in selecting the energy conversion and storage devices making up the system. The obtained LMES should be able to satisfy the fluctuating energy demand at all time and to minimize the total construction and operation cost of the system over its lifetime usually spanning several decades.

Many off-the-shelf numerical decision-aid tools already exist to assist people in the design of LMESs. However, most of these tools rely on strong assumptions and simplifications. For example, they assume that the capacity of an energy conversion device may take any value within a predefined continuous range, while this capacity should in fact be selected within a discrete list of values corresponding to the available models produced by equipment manufacturers. Moreover, these tools usually strongly limit the size of the considered instances to enable the resolution process to terminate within an acceptable computation time.

This PhD thesis focuses on the problem of optimally designing an LMES involving both energy conversion and storage devices. In terms of problem modelling, we improve the current state-of-the art in several directions. First, we allow to choose the capacity of the installed devices within a predefined discrete list of value. Second, we consider the fact that building an LMES is not a one-step but rather a multi-step process in which investment decisions are made little by little to adjust the system layout to the long-term increase of the energy demand. We thus seek to build a multi-phase strategic deployment plan for the LMES. Third, we incorporate in the objective function the operation cost of the system over its lifetime. To estimate this cost as accurately as possible, we build detailed daily operation schedules using hourly time steps for a set of representative days. These schedules take into account several realistic complicating features such as the partial load efficiency and the minimum working load of energy conversion devices.

The resulting optimization problem is formulated as a very large mixed-integer linear program. To solve it efficiently, we develop two new decomposition algorithms exploiting the specific bi-level structure of the problem. The first algorithm extends a previously published hierarchical decomposition algorithm, the second one is a generalized Benders' decomposition algorithm.

The proposed modelling and solving approach are applied to three real-life case studies located in China. Our numerical results first show that the proposed decomposition algorithms significantly outperform both the generic branch-and-cut algorithm embedded in a mathematical programming solver and the original hierarchical decomposition algorithm at solving the mixed-integer linear program to vi optimality. Our results also show that, even if some approximations are done in the problem modelling, the obtained deployment plans are of very good quality.

Résumé

Au cours des dernières décennies, l'industrie énergétique s'est efforcée d'améliorer l'efficacité de la production d'énergie, de réduire les émissions de gaz à effet de serre liées à la production et à la distribution d'énergie et de mieux intégrer les ressources énergétiques renouvelables. Les systèmes multi-énergies locaux (LMES) constituent une alternative intéressante pour atteindre ces objectifs ambitieux. Fondamentalement, un LMES est un système énergétique décentralisé produisant de l'énergie sous de multiples formes pour satisfaire les besoins énergétiques de clients situés dans son voisinage. Les clients correspondent à un ensemble de bâtiments appartenant, par exemple, à un campus, un complexe hospitalier ou un quartier urbain. Les LMES ont une efficacité de production plus élevée et un coût de maintenance plus faible que les systèmes énergétiques traditionnels fonctionnant au niveau d'un seul bâtiment.

Les LMES présentent donc de nombreux avantages pratiques. Toutefois, afin de fournir leurs meilleures performances potentielles, ils doivent être conçus avec soin. La conception d'un LMES consiste essentiellement à sélectionner les dispositifs de conversion et de stockage de l'énergie qui composent le système. Le LMES obtenu doit être capable de satisfaire à tout moment la demande énergétique fluctuante et de minimiser le coût total de construction et d'opération du système sur sa durée de vie, qui s'étend généralement sur plusieurs décennies.

Il existe déjà de nombreux outils numériques pour aider les gens à concevoir des LMES. Cependant, la plupart de ces outils reposent sur des hypothèses et des simplifications importantes. Par exemple, ils supposent que la capacité d'un dispositif de conversion d'énergie peut prendre des valeurs dans une gamme continue prédéfinie, alors que cette capacité devrait en fait être sélectionnée dans une liste discrète de valeurs correspondant aux modèles disponibles produits par les fabricants d'équipements. De plus, ces outils limitent généralement fortement la taille des instances considérées afin de pouvoir proposer des solutions dans un temps de calcul acceptable.

Cette thèse de doctorat se concentre sur le problème de la conception optimale d'un LMES impliquant à la fois des dispositifs de conversion et de stockage d'énergie. En termes de modélisation du problème, nous améliorons l'état de l'art actuel dans plusieurs directions. Premièrement, nous permettons de choisir la capacité des dispositifs installés dans une liste de valeurs discrètes prédéfinies. Deuxièmement, nous considérons le fait que la construction d'un LMES ne se fait pas en une seule étape mais est plutôt un processus en plusieurs étapes dans lequel les décisions d'investissement sont prises petit à petit pour ajuster le déploiement du système à l'augmentation à long terme de la demande d'énergie. Nous cherchons donc à construire un plan de déploiement stratégique en plusieurs phases pour le LMES. Troisièmement, nous incorporons dans la fonction objectif le coût d'opération du système sur sa durée de vie. Pour estimer ce coût aussi précisément que possible, nous construisons des plannings d'opération journaliers détaillés aux pas de temps horaires pour un ensemble de jours représentatifs. Ces plannings tiennent compte de plusieurs caractéristiques réalistes compliquées telles que l'efficacité à charge partielle et la charge minimale des dispositifs de conversion d'énergie.

Le problème d'optimisation qui en résulte est formulé comme un très grand programme linéaire à nombres entiers mixtes. Pour le résoudre efficacement, nous viii développons deux nouveaux algorithmes de décomposition qui exploitent la structure spécifique à deux niveaux du problème. Le premier algorithme étend un algorithme de décomposition hiérarchique précédemment publié, le second est un algorithme de décomposition de Benders généralisé.

L'approche de modélisation et de résolution proposée sont appliquées à trois cas d'étude réels situés en Chine. Nos résultats numériques montrent d'abord que les algorithmes de décomposition proposés sont plus performants que l'algorithme générique de branch-and-cut intégré à un solveur de programmation mathématique et que l'algorithme original de décomposition hiérarchique pour résoudre le programme linéaire en nombres entiers à l'optimalité. Nos résultats montrent également que, même si certaines approximations sont effectuées dans la modélisation du problème, les plans de déploiement obtenus sont de très bonne qualité.

Introduction

A Local Multi-Energy System (LMES), or Decentralized Energy System, produces energy within multiple forms in order to satisfy the energy needs of customers located in its vicinity. The produced energy can be in the form of heating power, cooling power, domestic hot water or locally generated electricity The served customers may correspond to e.g. a small set of buildings [Wak+19; Meh+12], a hospital complex [START_REF] Arcuri | A mixed integer programming model for optimal design of trigeneration in a hospital complex[END_REF] or a urban district [START_REF] Jennings | Modelling and optimization of retrofitting residential energy systems at the urban scale[END_REF].

LMESs contrast with individual energy systems in which each single building is equipped with its own heating or cooling system. LMESs are usually much more efficient than individual energy systems with respect to both energy conversion and cost. Namely, individual energy systems mostly rely on small-capacity energy production devices whose conversion efficiency is much lower than the one of the medium-capacity devices that can be used in an LMES. In addition, the energy consumption pattern of a building depends on its function. For instance, the energy consumption of an office building concentrates in weekday daytime while the one of a residential building is higher at night and during the weekend. An LMES serves buildings with various functions. As a consequence, the total demand that it will have to satisfy will be somewhat less variable than the one of an individual building and is less likely to fall to zero, which will reduce the idle time of the system. Finally, the local generation of electricity may contribute in reducing the dependence of the customers from the national and global electricity markets, in which the price of electricity displays strong fluctuations.

An LMES is thus a decentralized energy system involving a variety of energy conversion devices (photo-voltaic panels, combined heat and power generators, boilers, chillers. . . ) and in some cases, energy storage systems (hot/chilled water tanks, ice storage, batteries, . . . ). These devices are used to transform energy resources into multiple forms of energy to satisfy the demands of the system's customers. The energy resources may be local energy resources (solar, wind, waste heat, geothermal, . . . ) or conventional energy sources (electricity, gas, biofuels, . . . ). Figure 1.1 describes the structure of an illustrative local multi-energy system aiming at satisfying the local demand of power, gas, cold and heat by consuming and converting electricity and gas bought from the national distribution network.

Although LMESs display many practical advantages, they should be carefully designed in order to provide their best potential performance. This PhD thesis focuses on the LMES optimal design problem which can be described as follows. Given the available energy resources, the predicted total energy demand in the neighbourhood over the LMES lifetime and the available energy conversion and storage technologies, select the proper technologies to be installed in the system and FIGURE 1.1: Example of a local multi-energy system [START_REF] Xie | The optimal planning of smart multi-energy systems incorporating transportation, natural gas and active distribution networks[END_REF] determine the size (i.e. capacity) of the corresponding conversion and storage devices so as to minimize the total construction and operation cost of the system while satisfying the customers' demand at all time.

Indeed, design decisions will have a strong impact on the daily operational management of the LMES so that the system design should consider not only the construction cost but also the operation cost over the system lifetime which usually spans several decades. Namely, the operation cost typically represents between 20% and 60% of the total cost of LMES over its lifetime. Moreover, a given conversion device may be rather cheap in terms of purchasing and installation costs but lead to high operation costs due to e.g. a poor energy conversion efficiency.

The operation cost of an LMES, which corresponds mainly to the cost of the energy resources bought from external suppliers to power the installed conversion devices, is however challenging to compute. Namely, the energy demands are highly variable and display a daily, weekly and yearly seasonality together with random variations. The price of the energy resources to be bought is also time-varying. Moreover, due to the technical features of some of the conversion devices, the energy resource consumption of the system is not proportional to the produced energy. This implies that using an average or aggregate representation of the demand to estimate the operation cost may lead to a significant underestimation of its actual value. Thus, in order to accurately estimate it, a detailed schedule describing, on a hourly basis, the on/off status and the load allocation of each conversion device should be built for an horizon spanning a whole year.

Furthermore, an LMES is often built together with the district it will serve. This means that the forecasted yearly energy demands feature an upward trend over the years as the district is developed and new buildings progressively connect to the LMES. As a consequence, the design of an LMES should not be seen as a one-time decision but rather as a process in which investment decisions are made step by step, following the development of the district and the increase of the demand. This implies that a multi-phase strategic deployment plan should be built.

We thus need to simultaneously determine, for each investment phase of the deployment plan, the optimal system layout and the corresponding detailed operation plan ensuring that the customers' energy demands are satisfied at all time and computing an estimation of the operation cost. Therefore, the optimal design problem is also referred to as the simultaneous system design and operation optimization problem.

Context 1.2.1 Industrial context

This PhD thesis results from a collaboration between EDF China, EDF R&D France and the University Paris Saclay.

In addition to its activities related to the large-scale generation of electricity in nuclear reactors or coal-fired power plants, EDF China is involved in several projects, sometimes together with Chinese external clients and partner companies, to build and operate LMESs in China. For instance, EDF China have been operating since 2016, together with the Chinese company Datang Group, a district heating network in the city of Sanmenxia [EDF22a]. The operation of the network is based on the heat recovered from two existing Datang power plants and 2x17 km backbone pipes have been built to connect the thermal power plants to the distribution networks in three industrial and business districts of the city. The heating area represented a total of 12 millions m 2 in 2020. Similarly, EDF China is currently developing, in partnership with ChangFeng Energy, a multi-energy system in the city of Sanya [EDF22b]. This system will supply chilled water for air conditioning and will contribute to the supply of sanitary hot water to a set of buildings (hotels, shopping centers and hospital) located along the Haitang Bay, a fast-developing tourist area. It is estimated that the LMES will reduce CO2 emission by 20% as compared to using individual and isolated systems, which represents a reduction of 70,000 tons of CO2 emission each year.

These LMES projects usually start with a call for tenders launched by cities wishing to develop an LMES in one of their districts. Companies such as EDF China and its competitors may reply by presenting a tender describing potential layouts for the system, together with their costs. This tender process usually comprises two steps. The first step involves a global selection of the technologies to be used in the system but usually relies on a gross estimate of the system operation cost assuming e.g. that the operating cost of each device is strictly proportional to the amount of output power. The second step consists in providing a detailed description of the future system comprising the number and capacity of each type of conversion and storage device to be installed in this system. In this second step, an estimation, as accurate as possible, of the future operation cost of the system is required. Overall, the goal of EDF China while presenting a tender is to select and size the devices installed in the system in such a way to simultaneously ensure that all the energy demands will be satisfied and that the total design and operation cost will be as low as possible.

Currently, the detailed tender to be prepared for the second step of this process is built manually by EDF China practitioners. Basically, the cost of each considered potential system layout is estimated through a set of computations carried out in Excel files. These computations rely on strong assumptions on the operation of the system, (e.g. devices are activated to produce energy in a predefined order and and energy storing and releasing is limited to some specific time periods). The system layout to be presented in the tender is then selected through a manual trial-and-error process but there is no systematic exploration of the set of all potential layouts.

Clearly, there is a need for a decision-aid tool based on mathematical optimization to help EDF China practitioners identify the best system layout (together with the multi-phase deployment plan) to make sure that the LMES projects developed will be profitable.

EDF R&D has a strong experience in developing mathematical-programming based decision-aid to help in the strategic and daily management of energy systems. For instance, EDF Lab Chatou developed a software called Clevery to help in the daily operation of smart micro-grids and local energy systems. This system is able to build optimal operation schedules for a large set of small-scale systems whose layout is already known. It however does not allow to make design decisions for such systems.

Thus, EDF currently has no decision-aid tool dedicated to the long-term design of LMESs. Our work in this PhD thesis aims at developing such a tool.

Academic context

As discussed in the previous subsection, optimally designing an LMES is a very challenging problem for practitioners which may be overwhelmed by the exponential number of design possibilities to be considered and compared. The difficulty experienced by practitioners is supported theoretically by the fact that the LMES design problem was recently proved to be NP-hard in [GCW19], even when there is a single type of conversion device, no storage device and a single time-step to represent the demand at the operational level.

[KV21] provided a recent state of the art on the modeling and solution approaches proposed to design LMES. They highlighted the fact that mixed-integer linear programming (MILP) was one of the most used optimization techniques, probably thanks to its high flexibility in terms of problem modeling and to the availability of mature off-the-shelf commercial MILP solvers. The literature review focused on MILP models and decomposition approaches for the LMES design problem, to be presented in Chapter 2, shows that the research is this field is currently very active and that much progress has been made over the last two decades.

However, there remains a significant gap between the current state of the art on the LMES design problem and the needs of our industrial partner.

Indeed, we first note that all the reviewed works consider that the system is built in a single construction phase while, in our practical problem, the LMES is to be built and extended progressively following a multi-phase deployment plan. Second, in most previously published works, the sizing of technologies is oversimplified. Namely, these works assume that, for each technology, there is a single piece of equipment to be installed and that the capacity of this device may be set to any value within a continuous range. However, in practice, we may install several devices corresponding to the same technology, the capacity of which should be selected among a set of discrete values corresponding to the models available in the catalog of the equipment manufacturer. Third, due to the high complexity of the problem, the assessment of the operation cost is usually carried out by introducing some simplifications on the functioning of the conversion devices, in particular with respect to their conversion efficiency, and by using a limited number of representative days.

In terms of solution approaches, most of the reviewed papers rely on a direct resolution of the problem by an MILP solver. This direct resolution is only possible at the expense of some strong simplifications on the actual problem. Moreover, decomposition approaches exploiting the natural hierarchy between the design and operation decisions in the problem have been proposed. But their direct application on the problem under study here is not straightforward.

Clearly, some research is needed to extend the previously published works and to incorporate into the problem modeling the above-mentioned important practical features. Taking these features into account leads to an increase in the MILP size and consequently makes it harder to solve. Part of the work carried out in this thesis will thus be devoted to the development of efficient solution approaches for this problem.

Background on mixed-integer linear programming

In this work, we heavily rely on mathematical programming, more specifically on mixed-integer linear programming (MILP), to tackle the problem of optimally designing an LMES. We thus provide in this section a short background on MILP to ease the reading of the manuscript.

Definition

A mathematical programming problem is an optimization problem taking the following form [START_REF] Boyd | Convex optimization[END_REF]:

minimize f 0 (x) subject to f i (x) ≤ b i , i = 1, • • • , m x ∈ R n
The vector x = (x 1 , • • • , x n ) represents the set of optimization or decision variables, the value of which is to be determined by solving the optimization problem. The function f 0 : R n → R is the objective function and represents the performance criterion used to evaluate and compare the various solutions of the problem. There is a set of m inequalities defining the feasible space, i.e. the subset of vectors x which are considered as acceptable (feasible) solutions of the optimization problem. These set of inequalities is defined through constraint functions

f i : R n → R, i = 1, • • • , m and right-hand side values b 1 , • • • , b m .
In case the objective function and all the constraint functions are linear functions of x, the problem simplifies into a linear programming (LP) problem in the form:

minimize ∑ n j=1 c j x j subject to ∑ n j=1 a i,j x j ≤ b i i = 1, • • • , m x ∈ R n
Here, all input parameters c j , a i,j and b i are real numbers whose value is assumed to be known.

An LP problem can also be written in the matrix form, using the convention that the ≤ sign compares the two vectors located on each side element-wise:

minimize c T x subject to Ax ≤ b x ∈ R n
Linear programming is a mathematical optimization tool which is widely applied to design energy systems. This may be explained by the fact that it enables to flexibly model a rather large range of design problems and that it leads to the formulation of mathematical programs which are relatively easy to solve. However, linear programming presents some drawbacks and limitations and may lead to recommending system designs which are sub-optimal or even not relevant in practice. In particular, the fact that the decision variables x may take any continuous value poses some modeling problems to represent e.g. decisions that should take integer values (such as the number of energy conversion units of a given type that should be installed in the system) or logical conditions (such as the fact that if a device is turned on, it should produce an amount of energy above its minimum production load). The optimal solution of an LP progam may involve components of x with a frational value and simply rounding this frational value up and down does not always provide a feasible, let alone a good, solution of the optimization problem (see e.g. [START_REF] Wolsey | Integer programming[END_REF]).

We thus sometimes have to add constraints to the formulation stating that some decision variables should be binary or integer in the solution. If all the variables are restricted to be binary, we obtain a Binary Linear Programming (BLP) problem. If all the variables are restricted to be integer, the problem is an Integer Linear Programming (ILP) problem. If the formulation involves both real valued variables and integer/binary variables, the problem is called a Mixed Integer Linear Programming (MILP) problem. The variables taking resp. binary, integer or fractional values are often referred to as resp. binary, integer and continuous variables. Let p be the number of integer variables in the MILP problem with 0 < p < n, the matrix form of an MILP problem takes the following form:

minimize c T x subject to Ax ≤ b x ∈ Z p × R n-p

Formulation of a mixed-integer linear program

The formulation of a real-life optimization problem as an MILP problem, i.e. the transformation of an informal textual description of the problem into a mathematical model, can be done through e.g. the procedure described in Chapter 1 of [START_REF] Pochet | Production planning by mixed integer programming[END_REF]. This modeling relies on the definition of the indices, data, variables, constraints and objective function of the MILP problem.

Indices Similar objects are grouped into object classes represented by mathematical indices, allowing one to use indexed notation for data, variables, and constraints. For instance, the index y (resp. ϕ) may be used to represent a year within the set of years (resp. a phase within the set of investment phases) of the LMES lifetime. Similarly, the index m may be used to indicate a type of energy conversion device within the set of all considered energy conversion devices.

Data Data refer to the "known values" in the optimization problem. For instance, the installation price of an energy conversion unit is considered as an input parameter in our case, the value of which is assumed to be known before solving the optimization problem.

Variables Variables are the "unknown values" of an MILP model and represent the decisions to be made to define a solution to the optimization problem. When formulating real-life problems, for the sake of clarity and readability, we usually do not use a single variable vector x to represent all the decisions to be made, but rather a set of multi-dimensional vectors referred to by different names. Thus, in our context, the integer decision variable SD ϕ,m indicates the number of conversion and storage units of type m to be installed at the beginning of investment phase ϕ whereas STO t,m represents the amount of energy stored in the devices of type m at the beginning of the scheduling time-step t. A solution of an optimization problem thus consists in an assignment of a (binary, integer or continuous) value to each decision variable involved in the problem modeling.

Constraints Constraints are mathematical relationships between the input data and the decision variables and take the form of inequalities or equalities. Recall that in order to obtain an MILP problem, these relationships should be linear expressions of the decision variables. A feasible solution of an MILP problem is a solution in which the values assigned to the decision variables satisfy all the constraints. A solution is said to be infeasible as soon as it violates one or several constraints. The set of all feasible solutions is called the feasible region or feasible space of the optimization problem. If the feasible region of a problem is empty, the optimization problem is said to be infeasible.

Objective function

The objective function provides a way to evaluate or compare feasible solutions and to select the best or optimal solution among the feasible ones. Both its optimization direction and its mathematical expression should be clarified.

The optimization direction can be either maximization or minimization: these directions are interchangeable by adding a negative sign to the objective function. As for the mathematical expression, in an MILP problem, it should be a linear expression of the decision variables. Finally, if there exists a sequence of feasible solutions in the feasible region whose objective values are unbounded, in other words if the objective function of the problem can take any arbitrarily 'good' value, the problem is referred to as an unbounded problem. In some books, the objective value of a minimization problem is defined as +∞ if it is unfeasible and as -∞ if it is unbounded (and vice versa for maximimization problems).

Branch-and-bound algorithm

One of the reasons explaining the widespread use of mixed-integer linear programming to tackle optimization problem is the successful implementation of the branchand-bound (B&B) algorithm in mathematical solvers. The B&B algorithm basically carries out an enumeration of all feasible solutions of an MILP, but uses the information provided by the linear programming relaxation of the original problem to cut out solutions proven to be non-optimal whenever it is possible. Note that using a B&B algorithm does not guarantee to find the optimal solution of any MILP in a reasonable time. However, as compared with other resolution strategies, it provides a relatively efficient solving routine for a large class of MILPs. It is in any case a fundamental cornerstone and a good starting point for solving MILPs. Before introducing the details of algorithm, we will first explain a trivial but important proposition for MILPs.

Proposition 1.3.1. Let IP 1 := min{c T x : x ∈ F 1 } and IP 2 := min{c T x : x ∈ F 2 } be two MILP problems. We denote by F 1 and F 2 the feasible regions of IP 1 and IP 2 . Let z * 1 and z * 2 denote respectively the optimal objective value of IP 1 and IP

2 . If F 1 ⊆ F 2 , then z * 1 ≥ z * 2 .
Proof. The proof of this proposition is straightforward. If the optimal solution of IP 1 exists and can be found at point x * 1 ∈ F 1 , then by inclusion, x * 1 is also in F 2 and the optimal objective value of P2 is at most z * 1 . Moreover, if we adopt the convention that the optimal value of an infeasible problem is +∞ and the one of unbounded problem -∞, this proposition still holds. If IP 2 is infeasible, then F 2 is empty, which implies that F 1 is also empty and IP 1 also infeasible. On the other hand, if IP 1 is unbounded, then there exists a sequence of {x n } n in F 1 whose objective value is unbounded, i.e., lim inf n→∞ c T x = -∞. By inclusion, this sequence is also in F 2 , hence the problem IP 2 is also unbounded.

We then consider an MILP problem of the form IP := min{c T x : x ∈ F}, where F = P 0 ∩ Z p × R n-p and P 0 = {x ∈ R n : Ax ≤ b}. P 0 is the convex set obtained when removing the integrality constraints on the variables of IP. Problem LP 0 := min{c T x : x ∈ P 0 } is a linear program called the linear programming (LP) relaxation of the original MILP problem.

The B&B algorithm starts with the resolution of LP 0 . A solution of the LP relaxation is said to be integer feasible if it is feasible for the original MILP problem. There are three possible cases for the optimal solution of LP 0 : it may be infeasible, non-integer feasible or integer feasible. First, if it is infeasible, from Proposition 1.3.1, we deduce that IP is also infeasible. Second, if it is integer feasible, we can deduce from Proposition 1.3.1 that the solution of LP 0 is the optimal solution of IP. Third, if it is non-integer feasible, i.e., if it contains non-integral values for some integer variables, by Proposition 1.3.1, the optimal solution of LP 0 provides a lower bound of the optimal solution of IP.

In the first two cases, the algorithm may stop, having either proved that the MILP is infeasible or found an optimal integer solution. Let us now focus on the third case. We select one of the integer variables having a non-integral value in the LP relaxation and denote by x 1 this integer variable and by x1 ̸ ∈ Z its current value. We create two new LP problems: LP 1 (resp. LP 2 ) is obtained by adding constraint x 1 ≥ ⌈ x1 ⌉ (resp. x 1 ≤ ⌊ x1 ⌋) to the LP relaxation. The LP relaxation together with these two generated LP problems form a binary tree, which is usually referred to as the branch-and-bound search tree, and the operation of generating two new LP problems is called branching. The root node of the tree, indexed by 0, corresponds to Problem LP 0 whereas its immediate children, Node 1 and Node 2, are associated to the two LP problems LP 1 and LP 2 created by branching. Let P 1 and P 2 denote the feasible region respectively of these two new LP problems. P 1 and P 2 are two disjoint subsets of the feasible region of LP 0 , i.e., P 1 ∪ P 2 ⊆ P 0 and P 1 ∩ P 2 = ∅. Moreover, the corresponding integer sets

F 1 = P 1 ∩ Z p × R n-p and F 2 = P 2 ∩ Z p × R n-p form a partition of the feasible set F of IP, i.e., F 1 ∪ F 2 = F and F 1 ∩ F 2 = ∅. So the optimal solution of IP, if it exists, is either in F 1 or in F 2 .
We thus continue searching for the optimal solution of IP by repeating the process described above. At each iteration, we select an active node, branch from it and solve the LP problems corresponding to its newly created child nodes. For each child node, we again have three cases. If the solution found is not feasible, we prune the node and will not consider branching from it in the latter iterations. If the solution found is integer feasible, it means that it is a feasible solution of IP. We refer to the best integer feasible solution of IP found so far as the incumbent: each time a new feasible integer solution is found, its value is compared with the one of the current incumbent and if this value is better, the incumbent is updated. If the child node returns a non-integer feasible solution, we add it to the list of active nodes and start a new iteration by selecting another active node.

The B&B tree thus progressively grows with each branching operation. However, as such, the B&B algorithm is equivalent to a trivial enumeration of all integer solutions in the feasible region. In order to be more computationally efficient, we can decide to prune nodes corresponding to non-integer feasible solutions by comparing the objective value of the associated LP problem with the one of the current incumbent. If the objective value of the associated LP problem is worst (i.e. larger or smaller depending on the optimization sense) than the one of the current incumbent, we can prune this node and stop branching from it. Namely, all integer feasible solutions of IP located in the feasible sub-region of IP defined by the branching decisions made up to this node, if they exist, cannot be better than the current incumbent and thus cannot be the optimal solution. This directly relies on Proposition 1.3.1.

Proposition 1.3.2. In the B&B search tree of a minimization (resp. maximization) problem, the objective value of a parent node is always smaller (resp. greater) than the objective value of its child nodes.

Proof. The proof is also straightforward, since the feasible region of a child node is always a subset of the one of its parent node and from the Proposition 1.3.1, we show that this proposition is true.

The B&B algorithm terminates when the set of active nodes is empty, and the incumbent is the optimal solution of the original IP problem.

Algorithm 1 provides a detailed description of the branch-and-bound algorithm than can be used to solve for a minimization problem.

Branch-and-cut algorithm

The computational efficiency of a branch-and-bound algorithm heavily depends on the quality of the bound computed at each node of the search tree. In a nutshell, the better this bound, the sooner a branch can be discarded and the less nodes have to be explored before the algorithm converges. In the context of mixed-integer linear programming, this bound is most often obtained by solving the linear programming relaxation of the problem. The main advantage of this is that linear programs can be very efficiently solved by the simplex algorithm. However, in many cases, the linear programming bound is of poor quality, which negatively impacts the performance of the branch-and-bound algorithm. In order to solve this issue, branch-and-bound algorithms can be coupled with a cutting-plane generation approach, giving rise to a branch-and-cut algorithm. Basically, a cutting-plane generation approach aims at improving the quality of the linear programming bound by adding a set of linear inequalities to the problem formulation. These inequalities are chosen so as to cut away non-integer solutions that would otherwise be solutions of the linear relaxation and to improve the value provided by the linear relaxation by restricting its feasible space. The problem of finding a cut separating a non-integer solution from the feasible space of the linear relaxation of the MILP is called the separation problem. The branch-and-cut algorithms embedded in MILP solvers use generic cuts, i.e. cuts that may apply to any MILP or at least to a wide range of MILPs, such as Gomory fractional cuts, clique cuts or mixed-integer rounding cuts.

Thanks among others to the use of branch-and-cut algorithms, MILP solvers have undergone tremendous progress over the last thirty years. For instance, [START_REF] Bixby | A Brief History of Linear and Mixed-Integer Programming Computation[END_REF] reported that the machine-independent computational performance of CPLEX solver improved by a factor of 29000 between the 1.2 version of the software released in 1991 and the 11.0 version released in 2007. Yet, despite this progress, there still are many MILPs whose direct resolution by a mathematical programming solver leads to prohibitive computation times. In particular, this arises when the MILP to be solved involves a large number of integer variables and constraints, as this is the case in the LMES optimal design problem.

Research objectives and main contributions

In this thesis, we aim to create a decision-aid tool based on mathematical optimization to help EDF practitioners at optimally designing an LMES. This tool should take as input data a list of available energy storage and conversion technologies, to be potentially installed in the LMES, together with a set of time series describing the predicted evolution of energy resources' price and availability and of the clients' energy demands during the expected lifetime of the system. The output should be a multi-phase deployment plan indicating the technologies to be used in the LMES together with the number and size of the conversion and storage devices of each technology to be installed in each construction phase. The obtained plan should ensure that the system will be able to satisfy the clients' energy demands at all time and that the total cost of the system, which includes its design and operation cost, is minimized.

The numerical tool to be developed should be flexible enough to design a broad class of LMESs and thus should take into account a wide range of energy resources, energy demands and energy conversion and storage technologies. Moreover, it should make it possible to size the selected technology in a practically relevant manner, i.e. by selecting the conversion and storage devices to be installed in the LMES among the discrete list of models available in the catalog of the equipment manufacturer rather than by simply fixing the capacity within a continuous range. Finally, as the operation cost makes up a large part of the total cost of an LMES, the tool should be able to accurately estimate this operation cost. A large number of typical and extreme days should thus be used in the mathematical optimization model to represent, as best as possible, the various conditions under which the system will be operated. Moreover, it will be necessary to take into account some practical constraints related to the conversion devices, in particular a non-zero minimum working load and non-linear performance curves. Finally, the tool should be able to provide deployment plans with a reasonable computational effort. Namely, even if the design of an LMES is a long-term strategic decision, the tool might not be adopted by its end-users if the computation time needed to find a deployment plan exceeds a few hours.

As can be seen from the description of the state of the art to be provided in Chapter 2, none of the available off-the-shelf decision-aid tools is able to meet all these needs. Moreover, to the best of our knowledge, no previously published academic work investigated a mathematical optimization-based approach for the LMES design problem while simultaneously considering all the above-mentioned practical features. Our contributions are thus threefold.

The first contribution pertains to the modeling of the LMES design problem as a mixed-integer linear program. First, we explicitly consider the fact that building an LMES is most often not accomplished in a single step but that it is rather a process involving multiple construction phases, each of which expanding over one or several years, to adjust to the long-term increase of the energy demand. To the best of our knowledge, this is the first time the LMES design problem is modeled as a multi-phase investment problem. Second, our model allows to choose the capacity of the installed devices within a discrete list of available models. We thus not only decide whether a technology will be used or not in the LMES but also determine the number and capacity of the devices belonging to this technology that should be installed. Finally, to estimate the operation cost of the system as accurately as possible, we build detailed daily operation schedules using hourly time steps. Each of these operation schedules consider several realistic complicating features such as the partial load efficiency and the minimum working load of energy conversion devices. Clearly, taking into account all these aspects of the problem leads to the formulation of a huge mixed-integer linear program which as such is computationally intractable. We thus propose a modeling approach aiming at reducing the size of the obtained mathematical program and at removing the non-linearities. This approach relies on the use of a limited number of investment phases to make design decisions and on the careful selection of a rather large number of disjoint representative days to estimate the operation costs. Moreover, in order to obtain a mixed-integer linear program, we use a piece-wise linear approximation of the energy conversion performance curves in case these ones are non-linear. Finally, we propose to exploit the convexity of these performance curves to build aggregate operation schedules, i.e. schedules describing the number of active devices, the energy consumption and the energy production of each set of identical devices, rather than detailed operation schedules describing the on/off status, the energy consumption and the energy production of each individual device. This enables us to reduce to a significant extent the number of decision variables needed in the MILP formulation to build these operation schedules.

The second contribution consists in the application and extension of two decomposition approaches to efficiently solve the large-scale MILP problem. Both approaches exploit the natural hierarchy between the design decision variables (used to build the deployment plan) and the operation decision variables (used to build the operation schedules). This hierarchy translates into a special structure of the mathematical problem. Namely, once the design decisions determining the system layout at each investment phase are fixed, the problem decomposes into a set of independent operation sub-problems, one for each considered representative day, aiming at building an optimal operation schedule for this day under the current system layout. This special structure is exploited by the hierarchical decomposition method previously published in [START_REF] Yokoyama | Optimization of energy supply systems by MILP branch and bound method in consideration of hierarchical relationship between design and operation[END_REF] but our numerical results show that a direct application of this method on our multi-phase problem is not computationally efficient. We thus extend the algorithm of [START_REF] Yokoyama | Optimization of energy supply systems by MILP branch and bound method in consideration of hierarchical relationship between design and operation[END_REF] to improve its efficiency at solving our multi-phase problem. This is achieved in particular through the addition of single-phase no-good cuts into the upper-level problem to forbid deployment plans using a system layout found to be infeasible for a certain investment phase, and through the recording of the solution value of previously solved operation sub-problems to avoid repetitive calculations. Furthermore, we develop a new decomposition approach for the problem which can be seen as a generalized Benders' decomposition algorithm. Note that the classical Benders' decomposition algorithm requires that the second-stage sub-problems display a strong duality property, which is not the case here since the operation scheduling sub-problems include discrete variables. We thus customize the framework presented in [START_REF] Bolusani | A Framework for Generalized Benders' Decomposition and Its Application to Multilevel Optimization[END_REF] to extend the classical Benders' decomposition algorithm to a larger class of optimization problems. The proposed generalized Benders' decomposition algorithm relies on a new set of feasibility and optimality cuts. These cuts exploit the specific structure of the constraints coupling the first-stage and the second-stage decisions in the operation scheduling sub-problems and enable us to obtain an algorithm with a theoretically guaranteed finite and optimal convergence.

The third contribution deals with the evaluation of the proposed modeling and solution approaches on instances based on three case studies corresponding to three real-life LMESs currently under construction by EDF in China. Two of these LMESs are district cooling systems (DCS) aiming at meeting a demand for cold. The third one is a district cooling and heating system aiming at meeting a demand for cold and heat: for this system, the installation of a trigeneration system, i.e. of a device capable of converting gas into electricity, cold and/or heat, is considered. For each of these three projects, the deployment plan includes multiple construction phases. Moreover, each system potentially involves conversion devices with non-linear performance curves and a non-zero minimum load rate. Our numerical results showed that the two proposed decomposition approaches significantly outperform, in terms of solution quality and computation time, both an MILP solver directly solving the MILP formulation and the previously published hierarchical decomposition method at solving the problem. In particular, the proposed generalized Benders' decomposition algorithm was able to solve all the considered instances to optimality within two hours of computation. However, solving the MILP problem to optimality does not guarantee that the obtained deployment plan is optimal with respect to the initial optimization problem as some approximations were made during the problem modeling. We thus carried out a post-optimization simulation study to evaluate the quality of the deployment plans obtained for the two DCS projects. Overall, this study indicates that, provided the number of representative days used in the model is large enough, the design decisions do not depend on the subset of selected representative days and the operation cost is estimated with a good level of accuracy by the optimization model so that the obtained deployment plans may be confidently recommended for a practical implementation. This remainder of this manuscript is organized as follows. In Chapter 2, we review the academic works closely related to ours, focusing in particular on the ones using an MILP model and/or a decomposition-based solution approach. In Chapter 3, we provide a detailed description of the optimization problem. We start by introducing the main components of an LMES: the commodities and the energy conversion and storage technologies. We then explain a number of additional aspects to be considered when designing an LMES and formally state the resulting optimization problem in terms of input data, decisions, constraints and objective function. Chapter 4 is devoted to the problem modeling and to its mathematical formulation as an MILP. We present among others the modeling choices and assumptions we made to obtain an MILP of tractable size. We then provide a detailed description of the three case studies used to create instances of the problem. We finally present the numerical results obtained when attempting to directly solve these instances with CPLEX 12.8 solver using the proposed MILP formulation. In Chapter 5, we extend the hierarchical decomposition approach previously published in [START_REF] Yokoyama | Optimization of energy supply systems by MILP branch and bound method in consideration of hierarchical relationship between design and operation[END_REF] to make it more efficient at solving the large-scale MILP model formulated in Chapter 4. This approach is applied on the three case studies. Our computational results show a clear improvement of the computational efficiency on the instances related to the two DCS projects. However, out-of-memory issues were still encountered while trying to solve the instances related to the trigeneration project. In Chapter 6, a new generalized Benders' decomposition algorithm is proposed. This new generalized Benders' decomposition algorithm relies on the specific structure of the constraints coupling the first-stage and the second-stage decisions in the operation scheduling sub-problems and extends the scope of the traditional Benders' decomposition to the case where sub-problems are MILP problems. The results of our computational experiments, carried out on instances based on our three case studies, are reported. They show that the generalized Benders' decomposition algorithm clearly outperforms all the other considered algorithms at solving the problem. Moreover, the outcome of a post-optimization simulation study carried out to evaluate the quality of the deployment plans obtained for the two DCS projects is discussed. Finally, in Chapter 7, we provide a general conclusion and discuss several potential research perspectives.

Chapter 2

State-of-the-art on the LMES design problem

Introduction

We provide in this chapter an overview of the works closely related to ours. We first discuss in Section 2.2 papers proposing a mathematical-programming based model to optimize the design of an LMES. We also present several already existing numerical MILP-based tools that may assist in managing and designing an LMES. We then focus in Section 2.3 on the literature dealing with solution algorithms relying on a decomposition of the main problem into a strategic design master problem and one or several operation scheduling sub-problems.

Mathematical programming models

Several modeling and optimization techniques have been investigated to tackle the problem of optimally designing an LMES: see e.g. [START_REF] Klemm | Modeling and optimization of multi-energy systems in mixed-use districts: A review of existing methods and approaches[END_REF] for a recent review on existing approaches. Among all these techniques, mixed-integer linear programming appears to be one of the most widely studied. This may be explained by the fact that MILP offers a high flexibility for modeling the problem and that several off-the-shelf mature commercial solvers (e.g. CPLEX [START_REF]Cplex Optimizer[END_REF], Gurobi [Gur21] and Xpress [FIC21]) capable of solving large size MILPs to optimality with a reasonable computational effort are now available.

This section thus aims at presenting the current state of the art on MILP models for the LMES design problem. Our work initially focused on the optimal design of district cooling systems (DCS) supplying a single type of energy (cold) and was later generalized to a broader class of local multi-energy systems supplying energy under different forms (cold, heat, electricity...). In what follows, we thus first review the works related to the design of a DCS before broadening the discussion to papers investigating the design of more general LMESs.

Optimal design of district cooling systems

DCSs have been the subject of many research works. We refer the reader to the two recent literature reviews provided by [EA19] and [START_REF] Gang | District cooling systems: Technology integration, system optimization, challenges and opportunities for applications[END_REF] for a general introduction on this domain and concentrate on papers dealing with a mathematicalprogramming based approach to optimally design a central chiller plant.

As mentioned among others by [EA19], most previously published works focus on optimizing the DCS distribution network configuration. Thus, we could find only a limited number of papers discussing the design of the central chiller plant. [START_REF] Söderman | Optimisation of structure and operation of district cooling networks in urban regions[END_REF] studies the optimal design of a district cooling network in an urban region. His problem consists in simultaneously locating and sizing a set of electric-powered compression chillers and of cold water storage sites and in selecting cooling pipeline routes to distribute the chilled water to the customers' buildings. The optimization problem seeks to minimize the total cost of building the system and operating it, i.e. of running the installed chillers and pumping the water in the pipes. A set of 8 twelve-hour periods are used to represent the seasonal and daily (night/day) variations of the cooling demand. [KH15] considers a similar problem and uses an oriented graph to represent the potential DCS. This graph comprises a single supply node representing the central chiller plant whose location is assumed to be known, a set of demand nodes representing the customers' buildings and a set of connection nodes corresponding to pipe junctions. Their design problem involves choosing the chiller plant capacity and the storage tank capacity among a discrete set of options, in selecting the arcs of the networks where distribution pipes should be built and in sizing these pipes. [START_REF] Al-Noaimi | Optimal design of a district cooling grid: structure, technology integration, and operation[END_REF] extends the work of [START_REF] Khir | Optimization models for a singleplant District Cooling System[END_REF] to the case where chillers may be installed at different potential supply nodes in the graph. Both works assume that the cooling demand of each connected building is one-day periodic and stationary and use a single typical day divided into 4 to 24 scheduling periods to evaluate the operation cost of the system. [START_REF] Dana M Alghool | Optimization of design and operation of solar assisted district cooling systems[END_REF] addresses the problem of optimally designing a DCS based on an absorption chiller using the hot water obtained from a set of solar collectors to produce chilled water. The design decisions involve determining the type and area of the solar collectors, sizing the absorption chiller and the auxiliary boiler and choosing the capacity of the cold and hot thermal storage tanks. The operation cost of the system is estimated by considering 12 typical days (one per each month in the year) divided into 24 scheduling periods. In all the above-mentioned works, the obtained optimization problem is formulated as a mixed-integer linear program and directly solved by a mathematical programming solver.

Optimal design of local multi-energy systems

Unlike DCSs which supply energy under a single form (cold), an LMES involves multiple forms of energy and a wide variety of conversion and storage devices, which significantly complicates the design problem. Several articles ([Yok+15; WHY21; Zho+13b]) point out that directly solving to optimality the LMES design problem with an MILP solver might not be possible in practice. These numerical findings are backed by the recent theoretical results presented in [GCW19]. The optimal design problem, which is referred to as the synthesis problem in this paper, is proven to be NP-hard even if there is a single type of conversion device, no storage device and a single period to represent the demand at the operational level.

In what follows, we review papers dealing with the design of the central power plant of an LMES. Papers dealing with the structure of the pipeline networks, such as [START_REF] Carmen | An energy integrated, multi-microgrid, MILP approach for residential distributed energy system planning-A South Australian case-study[END_REF], [START_REF] Eugenia D Mehleri | Optimal design and operation of distributed energy systems: Application to Greek residential sector[END_REF] or [START_REF] Weber | Design and optimization of district energy systems[END_REF], are out of the scope of this study. Tables 2.1 and 2.2 provide a summarized description of a list of recent papers dealing with the optimal design problem of an LMES. Each work is described according to five main features. The first three features, which are related to the problem modeling, are presented in Table 2.1. We thus consider the structure of the LMES (i.e., the available energy resources, the form of the energy demands to be met, the available conversion and storage technologies), the temporal granularity (length of a scheduling time-step and number of representative days) used to estimate the operation cost and the representation of the conversion efficiency of the devices (constant or partial-load efficiency). The last two features, which are related to the mathematical formulation and the solution approach, are presented in Table 2.2. We thus report for each paper the type (continuous, integer or binary) of the design and operation variables, the characteristics of the mathematical formulation (mixed-integer linear or non-linear programming) and the type of solution approach (direct resolution by a solver, evolutionary algorithm, Benders' decomposition...). Notice that in [START_REF] Zhou | A two-stage stochastic programming model for the optimal design of distributed energy systems[END_REF] and [START_REF] Steen | Modeling of thermal storage systems in MILP distributed energy resource models[END_REF], the proposed mathematical model can be used to design a wide range of LMESs and does not assume a specific structure of the system: this is indicated by noting "All" for the set of energy resources, energy demands and conversion technologies in Table 2.1.

Tables 2.1 and 2.2 show among others the wide variety in terms of structure of the LMESs investigated in the literature. In what follows, we review the papers while focusing mainly on the way the structure of the LMES is determined in the mathematical model and on three aspects which have a strong impact on the formulation of the operation sub-problems: the representation of the performance curves, the incorporation of an imposed minimum load for conversion devices and the presence of storage devices. The solution approaches will be compared and discussed in the next section.

We first compare the papers with respect to the way design decisions about the type and capacity of the installed devices are made: see column 'Design variables' in Table 2.2. We note that most papers use a combination of binary and continuous ('B+C') variables. The binary variables are used to decide whether a specific technology (heat pump, electric heater, gas-fired boiler, photo-voltaic panels...) or the device corresponding to this technology is selected or not. Continuous variables are used to determine, for each selected technology, the capacity of the installed device. Note that some papers [Zho+13b; UDT19; Ste+15; AFF07] use only continuous variables for the system design. However, sizing a technology, i.e. determining the capacity of the corresponding device installed in the system, using continuous variables is not relevant for many industrial case studies. Namely, in practice, the capacity of the installed devices cannot be set to any arbitrary value but rather has to be chosen within a set of discrete values. This set is defined by the catalog of available models provided by the equipment manufacturer. Moreover, using binary variables to decide whether a technology will be used or not in the LMES relies on the assumption that there is a single type of device available for this technology whereas there may be several types of devices, each one corresponding to a given capacity, for a given technology. This is why integer variables indicating the number of devices of each type (i.e. of each technology and each possible value of the capacity) to be installed should be used to determine the structure of an LMES. We found only four papers [Zho+13a; Ren+10; CSL11; Yok+15] using integer variables to determine the number of conversion devices to be installed in the system. Moreover, in [START_REF] Omu | Mixed integer linear programming for the design of solar thermal energy systems with short-term storage[END_REF], a single integer variable is used to determine the number of solar thermal collector panels to be built, while the size of the electric heater is chosen within a continuous range.

The second important aspect when modeling the LMES design problem by MILP is the representation of the performance curves of the conversion devices: see column 'Conversion efficiency' in Table 2.1. Recall that these curves provide the amount of energy consumed by a device as a function of the energy produced. The conversion efficiency of a device may depend on many factors such as the ambient temperature (see e.g. the solar thermal collectors in [START_REF] Omu | Mixed integer linear programming for the design of solar thermal energy systems with short-term storage[END_REF]) or the capacity of the device (see e.g. the gas turbine of the combined heat and power unit in [START_REF] Weber | Design and optimization of district energy systems[END_REF]). As long as these factors do not depend on the scheduling decisions, the conversion efficiency can be considered as constant, which results in linear performance curves.

However, for some types of device (e.g. electric chillers), the conversion efficiency depends on the load rate of the device, which is a scheduling decision, and most often, its partial-load efficiency is smaller than its full-load or nominal efficiency. This translates into non-linear performance curves, which are difficult to handle in the MILP formulation. Note that the vast majority of the reviewed papers assume a constant conversion efficiency and use linear performance curves. We could find only three exceptions. In [START_REF] Zhou | Impacts of equipment off-design characteristics on the optimal design and operation of combined cooling, heating and power systems[END_REF], the non-linear performance curves of the gas boiler, the gas engine and the absorption chiller are linearized using a piece-wise linear approximation approach. The fuel consumption of the combined heat and power (CHP) component in [START_REF] Elsido | Twostage MINLP algorithm for the optimal synthesis and design of networks of CHP units[END_REF] depends on both the load rate and the extraction valve opening rate: the authors use a triangle method [START_REF] Ambrosio | Piecewise linear approximation of functions of two variables in MILP models[END_REF] to approximate this function by a two-dimensional piece-wise linear function. In [START_REF] Urbanucci | A comprehensive methodology for the integrated optimal sizing and operation of cogeneration systems with thermal energy storage[END_REF], the partial-load efficiency of a combined heat and power unit is explicitly considered in the model, leading to the formulation of an MINLP.

There is a third feature strongly impacting the problem modeling: the minimum load rate of the conversion devices. Namely, in some cases, the output power of a device cannot take any arbitrary value between zero and its maximum capacity. Instead, a device, when active, has to produce a minimum amount of power. As can be seen from the column 'Operation variables' of Table 2.2, most papers use binary variables to represent the on/off status in each scheduling time-step of the single piece of equipment installed for each technology. However, models using a binary variable for selecting a technology and a continuous variable for sizing it neglect the fact that there will in fact be multiple devices of this technology in the system and can therefore not optimize the detailed operation of these multiple devices. The simulation study presented in [START_REF] Chen | Design and operation optimization of multi-chiller plants based on energy performance simulation[END_REF] shows that the actual energy consumption of a DCS designed while explicitly building optimal schedules for each of the chillers involved in the system can be around 20% that the one of a DCS designed while using a simplified control strategy not taking into account the presence of multiple chillers in the system. Among the three papers using integer design variables to determine the number of devices of each type to be installed, we may distinguish two cases.

[Zho+13a] builds a detailed operation schedule for each individual device thanks to the use of binary operation variables whereas [START_REF] Hongbo Ren | Multi-objective optimization for the operation of distributed energy systems considering economic and environmental aspects[END_REF] builds an aggregate operation schedule for the set of similar devices thanks to the use of integer operation variables. As for [START_REF] Carvalho | Optimal synthesis of trigeneration systems subject to environmental constraints[END_REF], the authors' focus is mostly on the environmental impact of the system design and operation details, such as the minimum load rate of the conversion devices, are not considered in their model.

Finally, the last aspect to be discussed is the presence of energy storage devices in the LMES: see column 'Storage technology' in Table 2.1. When there is no storage device in the system, the only way to satisfy the energy demand in a given time-step is to produce it through a conversion device during this time-step. Therefore, once the system design is determined, each operation sub-problem consists in ensuring the demand satisfaction in a single time-step and is thus of limited size. In contrast, the presence of an energy storage device results in a coupling between the scheduling time-steps as the amount of energy stored in the device at the end of a time-step defines the amount of energy available in this device at the beginning of the following time-step. These means that the single time-step operation sub-problems mentioned above are linked together by inventory balance equations, leading to larger and more difficult to solve operation sub-problems. As shown in Table 2.1, about one third of the reviewed articles do not involve energy storage in their system. Among the articles considering energy storage, the most common method to assess the operation cost of the LMES consists in using representative days, each one typically covering 24 one-hour time-steps. These representative days are chosen to represent as best as possible the energy demand profile along the year and are assigned a weight corresponding to the number of actual days in the year they represent. [START_REF] Morvaj | Optimization framework for distributed energy systems with integrated electrical grid constraints[END_REF] and [START_REF] Wakui | Two-stage design optimization based on artificial immune system and mixed-integer linear programming for energy supply networks[END_REF] both use a set of disjoint representative days, i.e. a set of representative days not linked to one another by inventory balance equations, and build a separate schedule for each of these days. Note that they consider that the entering energy inventory at the beginning of a representative day should be equal to the leaving energy inventory of the same day. This amounts to assuming that each representative day will be cyclically repeated during the year a number of times equal to its weight. [UDT19] consider inventory balance equations for all time-steps of the year and apply the rolling-horizon approach described in [START_REF] Bischi | A rolling-horizon optimization algorithm for the long term operational scheduling of cogeneration systems[END_REF]. This approach solves a sequence of operation sub-problems spanning two consecutive days. At iteration k, an operation schedule is buit for two consecutive days, denoted by Day k and Day k + 1 and the inventory level at the end of Day k is recorded. This level is used during iteration k + 1 to define the inventory level at the beginning of the Day k + 1 and build an operation schedule for Days k + 1 and k + 2. [Gab+18] also build an operation schedule spanning a whole year and taking into account all the inventory balance equations. However, they use an approximation enabling them to significantly reduce the number of binary scheduling variables and to obtain a computationally tractable problem. Namely, the authors propose to first build clusters of days. Then, when building the operation schedule for the whole year, they impose that the binary scheduling variables describing the on/off status of each device in each time-step should take the same value for all the days belonging to the same cluster. They however allow the continuous variables describing e.g. the output power of each device to vary from one day to the next in each cluster.

We end this section by a short discussion on the available decision-aid tools that may be used to optimize the cost of energy systems ([VBGS15; Man14; Lop+18]). Noe of them seems to be versatile enough to model and optimize all the aspects of an LMES that we wish to take into account. For example, the tools Balmorel [Bal21; Wie+18] and EnergyPLAN [Ene21]) can only optimize the operational cost. The tool DER-CAM ([DC21]) models an LMES by an MILP and can find an optimal strategy for the planning and operation of the system. It however does not allow to make design decisions and cannot take into account the fact that some devices may be operated under multiple modes (such as a dual-mode electric chiller producing either cold water or ice). Finally, eTransport ([BSW07]), another tool based on MILP, focuses on the topology structure of the system and the timing of the investments to be made. But its modeling of the operational functioning is not detailed enough to be used in our case as it considers neither the partial load efficiency, nor the minimum load rate, nor the multiple operation modes of energy conversion devices.

Decomposition-based solution approaches

Modeling the optimal design of local or district energy systems as a combinatorial optimization problem leads to the formulation of large-size mathematical programs which are very challenging to solve. Consequently, a wide variety of operation research-based solution approaches have been investigated for these problems. We focus here on the subset of approaches exploiting the natural hierarchy between the design and the operation decisions to decompose the optimization problem. More precisely, in the considered works, the optimization problem is split into a master problem determining the general structure of the energy system and one or several 2.2, as well as [START_REF] Fazlollahi | Multiobjectives, multi-period optimization of district energy systems: II-Daily thermal storage[END_REF] which is not in the table. These works consider a master problem in which all the operational decision variables and constraints have been removed and explore the set of feasible design decisions through a scatter search algorithm [UDT19; Zho+13b; Els+17] or a multi-objective evolutionary algorithm [START_REF] Weber | Design and optimization of district energy systems[END_REF][START_REF] Fazlollahi | Multiobjectives, multi-period optimization of district energy systems: II-Daily thermal storage[END_REF]. The operation cost of each potential design solution considered during this exploration is evaluated by solving a set of scheduling sub-problems using a mixed-integer linear programming solver.

Benders' decomposition algorithms such as the ones investigated e.g. by [START_REF] Falke | Investment planning of residential energy supply systems using dualdynamic programming[END_REF] and [START_REF] Li | Optimal design and operation of energy systems under uncertainty[END_REF] also consider master problems comprising only design decisions. The impact of these decisions on the operation scheduling sub-problems is taken into account in the master problem through a set of feasibility and optimality cuts. Note however that the generation of these cuts rely on the assumption that the scheduling problems are linear or convex programs involving neither binary nor integer variables. This implies that these algorithms may not straightforwardly extended to incorporate features such as a minimum working load or multiple operating modes for the conversion devices. [WHY21] also proposes to apply Benders' decomposition to a design problem in which sub-problems contains binary variables. The feasibility cuts are added in the form of no-good cuts but the optimality cuts are generated by using the dual solution of the linear relaxation of the scheduling sub-problems. This means that the operation cost is under-estimated in the master problem and that the obtained design solution is sub-optimal.

A decomposition approach based on the selection of representative days is proposed in [START_REF] William W Tso | A hierarchical clustering decomposition algorithm for optimizing renewable power systems with storage[END_REF]. The system design is determined by solving a MILP problem with a small number of representative days, and then, its corresponding operation cost over the full time horizon is computed. If the gap between the total cost computed with representative days and the total cost over the full time horizon is larger than the threshold, the algorithm increases the number of representative days and repeats the whole process.

Finally, several works consider master problems in which some or all of the operational decision variables and constraints are kept within the master problem but only (relaxed) continuous operation variables are considered. In this way, [START_REF] Iyer | Synthesis and operational planning of utility systems for multiperiod operation[END_REF] proposes a bi-level decomposition method. At the upper level, they formulate a master problem corresponding to the initial optimization problem in which all the binary operation variables have been removed but the continuous ones are kept. This master problem selects the conversion devices to be installed. Each time a potential system design is found, a set of lower level sub-problems is solved, each one taking as input the current system design and seeking to build optimal singleperiod schedules. The authors introduce the concept of design cuts to tighten the gap between the solutions obtained at both levels. [Yok+15] also presents a bi-level decomposition method. The upper level problem is a relaxed version of the initial optimization problem, in which all operation integer variables are kept but relaxed to be continuous. This problem is solved by a Branch & Cut algorithm. Each time a potential incumbent solution is found during the tree search, a set of single-period independent scheduling sub-problems is solved to check the feasibility and accurately compute the operation cost of the current system design.

Conclusion

We provided in this section an overview of the literature closely related to our problem. This overview enabled us to identify a significant gap between the current state of the art on the LMES design problem and our industrial problem.

Namely, we first note that all reviewed works consider a single construction phase while, in our practical problem, the LMES is to be built and extended progressively following a multi-phase deployment plan. Second, in most previously published works, the sizing of technologies is oversimplified. Namely, these works assume that, for each technology, there is a single piece of equipment to be installed and that the capacity of this device may be set to any value within a continuous range. As mentioned above, in practice, we may install several devices corresponding to the same technology, the capacity of which should be selected among a set of discrete values corresponding to the models available in the catalog of the equipment manufacturer. Third, due to the high complexity of the problem, the assessment of the operation cost is usually carried out by introducing some simplifications on the functioning of the conversion devices, in particular with respect to their conversion efficiency, and by using a limited number of representative days. In terms of solution approach, most of the reviewed papers solve the MILP problem directly by a mathematical solver. Some articles propose decomposition algorithms exploiting the hierarchical structure of the LMES but only a handful of these algorithms has a guaranteed finite and optimal convergence. This is why, in this thesis, we aim to extend the current state of the art in order to meet the needs of our industrial partner. We thus seek to develop an MILP model of the LMES design problem which can simultaneous take into account (i) multiple construction phases, (ii) a wide range of energy resources, energy demand, energy conversion and storage devices, (iii) a realistic choice of these devices among the models available at the manufacturer, (iv) the non-linearity of the performance curves of some conversion devices and (v) a number of representative days sufficiently large to obtain a good estimation of the actual operation cost of the system. In order to solve the resulting large-size combinatorial optimization problem to optimality within an acceptable computation time, we will investigate and extend decomposition algorithms exploiting the bi-level nature of the problem.

Chapter 3

Problem description

Introduction

This chapter provides a detailed description of the optimization problem under study and introduces the various elements to be taken into account when modeling it as a mathematical program.

In Section 3.2, we first introduce the key components of a generic LMES and provide a short technical description for the main energy conversion and storage technologies encountered in LMESs. We then further explain in Section 3.3 a number of additional aspects to be considered when designing an LMES. In Section 3.4, we formally state the optimization problem by defining the problem input data, the decisions to be made, the main constraints to be respected and the performance criterion to be optimized. Section 3.5 finally describes three application case studies corresponding to LMESs EDF is currently developing in China.

Components of an LMES

An LMES usually comprises one or several local energy plants, producing energy under multiple forms, and a distribution network to distribute this energy to the buildings connected to the system. In the present work, we will consider the case where the LMES involves a single local energy plant and will focus on choosing and sizing the conversion and storage devices in this energy plant. The problem of optimally designing the distribution network is thus outside of the scope of this PhD thesis. This corresponds to situations, often encountered in practice, in which the LMES system to be designed will use local pipeline networks which either already exist or will be built by a partner company. Such situations were considered e.g. in [AFF07; GCW19; Sch+18].

The central energy plant of an LMES involves a large number of various components and the complexity of its structure makes the design problem difficult to solve. However, not all the components have the same importance when making long-term decisions on the system layout. Our focus will be on choosing and sizing the energy conversion and storage devices so that the structure of the central energy plant can be significantly simplified without a big loss of accuracy. Thus, auxiliary components (such as pipes, valves, pumps and etc.) will not be considered in our problem modeling. Namely, the cost of these auxiliary elements is usually either fixed or proportional to the cost of the selected conversion and storage devices. It can thus be incorporated in the construction cost of these devices without major consequences on the long-term design decisions. Therefore, we will consider two main classes of components of the LMES: commodities and technologies. Commodities refer to the various forms of energy involved in the system: electricity, cooling power, heat... Technologies are the processes which transform one or several types of commodities into one or several other types of commodities.

In what follows, we discuss in more detail the commodities and technologies to be considered.

Commodities

We first classify the available commodities into three sub-types according their role in the LMES: resource commodities, supply commodities and intermediate commodities.

A resource commodity refers to a form of energy consumed by the system (e.g electricity from the national grid, natural gas and solar radiation). It is associated to a buying price, and/or to a limit on the available instantaneous power.

A supply commodity refers to a form of energy produced by the LMES. It can be used to satisfy the customers' demand (e.g. cooling power, heating power, electricity), in which case it is associated to a time-varying demand that should be satisfied at all time. In some cases, a supply commodity can also be sold (e.g. generated electricity sold to the national grid) and it is associated with a selling price.

Finally, an intermediate commodity refers to a form of energy produced and consumed within the system (e.g. hot water for storage, locally generated electricity, ice, exhausted gas). Intermediate commodities, as they exist only within the system, should be balanced between production and consumption.

Besides this classification based on their role in the LMES, commodities can also be described according to their nature. We introduce in what follows several types of commodities commonly used in LMESs.

Electricity An LMES generally consumes electricity coming either from the national grid or from a local generation. According to our classification, the electricity bought from the national grid is a resource commodity for which a time-varying price per unit of consumed electricity is charged. In some cases, there is also a limit (determined e.g. by a contract with the electricity supplier) on the electric instantaneous power that the LMES is allowed to consume. Electricity can also be generated locally by Combined Heat and Power (CHP) units or photo-voltaic (PV) panels. This locally generated electricity can be seen as an intermediate commodity consumed within the LMES by some conversion devices and/or as a supply commodity sold to the national grid or used to satisfy the customers' demand for electricity.

Natural gas Natural gas is a resource commodity for an LMES. The unit buying price of natural gas is usually much more stable than the one of electricity. The flame temperature of natural gas is around 1980 • C, which means that its combustion generates a lot of heat. Natural gas can thus be used to supply an internal combustion engine generating electricity and/or to provide direct heating. When it is used to generate electric power, the exhaust heat can be collected by heat exchangers or absorption chillers to produce heating or cooling power. This technique is widely used to improve the energy efficiency of natural gas combustion.

Chilled water

The thermal energy produced in the central plant of an LMES is transferred to the customers via a pipeline network, in which the circulating medium is usually water. Chilled water has a temperature lower than the one of the ambient environment and is used among others to provide air conditioning to the customers in summer. Chilled water is obtained either by cooling down water through devices such as chillers and heat pumps or by melting ice held in a thermal storage tank. Chilled water is transferred to the customers' buildings and is thus considered as a supply commodity with a given time-varying demand to be satisfied.

Ice Similar to chilled water, ice is produced to provide cooling power to the customers. Ice is obtained by cooling down water below its melting point through specific devices such as dual-mode chillers. Once produced, ice is usually stored for a short period of time (typically a few hours) in a thermal storage tank. When necessary, the LMES operator can melt the stored ice to obtain chilled water to be transferred to the customers' buildings. Ice is thus a way of storing thermal energy and is not used directly to satisfy the customers' demand. In an LMES, it is considered as an intermediate commodity.

Hot water

In winter, hot water is circulated in the pipeline network to satisfy the customers' heating demand. The heating energy carried by hot water can be obtained from the combustion of gas (e.g. in a CHP unit or a boiler) or from an electricity-driven thermodynamic cycle (e.g. in an air source heat pump). The hot water may be directly circulated to the customers' side to provide heating energy or be stored for a short period of time in a thermal storage tank. Hot water is thus both an intermediate and a supply commodity in an LMES.

Solar radiation Solar radiation is seen as a clean and sustainable source of energy but it has some practical drawbacks. First, it is of course only available at daytime. Second, its availability strongly depends on the weather conditions and on the season. In practice, researchers use local historical solar radiation records to forecast the future availability of solar energy: see e.g. [Alg+20; Tso+20]. Solar radiation is thus considered as a resource commodity which is free (the unit buying price is zero) but has a strongly time-varying limited availability. Solar radiation can be converted to electricity by photo-voltaic panels and/or to hot water by solar collectors.

Technologies

The energy conversion and storage technologies that may be used in an LMES differ among others with respect to their technical functioning and to the consumed and produced commodities. Due to the high variety of available technologies, it is not possible to provide here an exhaustive description. We thus only briefly present some technologies commonly used in LMES, focusing in particular on the ones encountered in our application case studies. For each considered technology, we present among others the consumed commodities, the produced commodities, the working load range, the conversion efficiency and the sizing pattern.

Electric chillers An electric chiller is a machine converting electricity into cooling power. This cooling power can be in the form of chilled water, which circulates through a pipeline network to absorb the heat in the buildings served by the LMES, or in the form of ice, which is stored in a thermal storage tank for a few hours and melted into chilled water afterwards. Not all chillers produce ice. We thus refer to chillers producing only chilled water as single-mode electric chillers (SMEC) and to chillers producing both chilled water and ice as dual-mode electric chillers (DMEC). An electric chiller has a working range defined by a maximum load, i.e. a maximum amount of cooling power it can produce per hour, and a minimum load, i.e. a minimum amount of cooling power it produces per hour when turned on. In practice, the minimum load under which it cannot be operated is equal to 10% to 20% of its maximum load. The efficiency of a chiller is usually characterized by its coefficient of performance (COP), defined as: COP = removed heat input electricity . The value of this COP is not constant: it varies with the chiller load, i.e. with the amount of produced cooling power, and with the ambient temperature. For example, at an ambient temperature of 30 • C, the COP of a single-mode chiller is about 6 when it is operated at full load, but goes down to 4 (resp. to 3) when the chiller is operated at a 30% (resp. 20%) load rate. Taking into account this varying part-load efficiency is especially important when the cooling demand is satisfied by multiple chillers running in parallel: the load rate of each active chiller should be carefully chosen to maximize the overall conversion efficiency. When designing an LMES in which some cooling capacity is required, we have to select the types of chillers to install by choosing among the available types of chillers provided in a given predefined list. Each type of chiller corresponds to a given category (either SMEC or DMEC), to a given working range (i.e. a minimum and maximum load) and to a given set of performance curves. These performance curves give the amount of consumed electricity as a function of the produced cooling power for some discrete values of the ambient temperature. The design decisions thus correspond to determining the number of chillers of each type to be installed in the system. CHP unit A Combined Heat and Power (CHP) unit, also called cogeneration system, converts natural gas to electricity and heating power. It contains two main components: a gas engine and a heat exchanger. A gas engine is an internal combustion engine which consumes gaseous fuel (usually natural gas in an LMES) to generate electric power. The exhaust heat is collected by the heat exchanger to obtain heating power. As mentioned above, the locally generated electricity is an intermediate commodity consumed by other devices in the LMES and/or a supply commodity distributed to the customers and sold to the national grid. The heating power produced is either an intermediate commodity stored in a thermal storage tank or a supply commodity used to satisfy the heating demand of the customers. The gas engine has a minimum working load which is usually set to more than half of its maximum working load or capacity [START_REF] Bischi | A detailed MILP optimization model for combined cooling, heat and power system operation planning[END_REF]: this means that a CHP unit should be operated within a minimum and maximum working load. Moreover, since the thermal energy is collected from the waste heat of power generation, heating power and electric power are always produced simultaneously in a CHP unit. A CHP unit is thus associated with both a heat efficiency coefficient and an electricity efficiency coefficient measuring the conversion efficiency of gas into electricity and heat. These two coefficients are usually constant, i.e. their value depends neither on the ambient temperature nor on the CHP load. When installing a CHP unit in an LMES, we have to select the type of CHP unit among a discrete list of available types given in a catalog, each type of CHP unit corresponding to a given gas engine and a given heat exchanger. The design decisions thus correspond to determining the number of CHP units of each type to be installed in the system.

CCHP unit A Combined Cooling, Heat and Power (CCHP) generation unit, also called trigeneration system, is a CHP unit linked with an absorption chiller. The absorption chiller uses some of the heat produced by the CHP unit to produce chilled water for air conditioning. A CCHP unit thus converts natural gas into electricity, heating power and cooling power. A CCHP unit cannot simultaneously produce heating and cooling power. It should be either on heating mode (i.e. producing heating power and electricity) or cooling mode (i.e. producing cooling power and electricity). Similar to a CHP unit, a CCHP unit is associated with a heat efficiency coefficient, a cooling efficiency coefficient and an electricity efficiency coefficient. All these three coefficients are constant as they vary neither with the ambient temperature nor with the load rate. Moreover, a CCHP unit must be operated within a minimum and maximum working load. When installing a CCHP unit in an LMES, we need to select the type of CCHP unit in a list of available types. Each type refers to a specific type of gas engine, a given heat exchanger and an a given type of absorption chiller, which are linked as a whole. The design decisions thus consist in determining the number of CCHP units for each type to be installed in the LMES.

Boiler A boiler converts natural gas to heating power. Basically, it is a pressure vessel in which water is heated through the combustion of natural gas. This hot water then circulates through the pipeline network to transfer heating power to the customers' buildings. Boilers may be operated at any value between 0 and their maximum working load or capacity, i.e. the minimum working load is equal to 0. Their conversion efficiency is constant and does not vary with the load or the ambient temperature. Similar to chillers and CHP units, installing a boiler in an LMES requires to select a type of boiler from a list of available models in a given catalog. The design decisions thus correspond to determining the number of boilers of each type to be installed in the system.

Heat pump A heat pump converts electricity to heating or cooling power. The thermal energy is transferred from the customers' buildings to the environment thanks to the circulation of a refrigerant which is driven by electricity. Heat pumps differ with respect to their environmental thermal source: there are e.g. air source heat pumps, geothermal heat pumps and water source heat pumps. In contrast to electric chillers which provide only cooling power, a heat pump can work either in heating mode or in cooling mode by reversing the circulation direction of the refrigerant. The cooling power produced by a heat pump cannot be stored and must be directly consumed by customers, while the heating power can be stored in the form of hot water. A heat pump may be operated at any value between 0 and its maximum working load or capacity. For a given (heating or cooling) mode, its conversion efficiency, which is defined as the ratio between the thermal energy production and the electricity consumption, varies only with the ambient temperature and does not depend on the load rate. When designing an LMES, we are given a discrete list of available heat pump sizes and we should choose the number of heat pumps of each size to install.

Thermal storage Thermal storage in an LMES usually refers to a container of heat or cooling material such as hot water or ice. There are many reasons to integrate thermal storage in an LMES. First, thermal storage may be useful to handle the shortterm (i.e. intra-day) variations of the energy demand. Thus, the system can produce more thermal energy than what is requested by the customers in low-demand periods, store the excess in the thermal storage container for a few hours and use this excess later in peak periods in which the customers' demand is high. Thus, a thermal storage may enable us to avoid installing energy conversion devices with a very large capacity to be used only in a small number of peak periods. Second, thermal storage may also be useful in an LMES involving co-generation devices such as a CHP unit. Namely, recall that with this technology, heat power and electric power are generated simultaneously when the gas engine is turned on. But the customers' thermal demand and electric demand do not necessarily vary in the same way. In this case, the surplus heat energy produced by the CHP unit can be stored in the thermal storage and used later instead of being rejected in the environment [START_REF] Bischi | A detailed MILP optimization model for combined cooling, heat and power system operation planning[END_REF]. A thermal storage may be used at any value between 0 and its maximum storage capacity. When designing an LMES, this storage capacity can usually be flexibly chosen to any value within a predefined continuous range. At the operational level, two main constraints should be complied with when operating a thermal storage: the amount of stored energy should not exceed the chosen thermal storage capacity and the evolution of the stored quantity should respect energy inventory balance equations stating that the energy flows into and out of the container together with the energy stored within the container respect the law of energy conservation.

RES diagram and superstructure

The design of an LMES starts from the high-level description of the potential energy system through a tool commonly used in the field of energy systems management called a Reference Energy System (RES) diagram. This diagram describes the commodities, the energy conversion and storage technologies, the energy flows between these technologies and the energy inputs and outputs. The commodities are represented by vertical segments and the technologies by rectangles. Arrows are used to indicate the commodity flows within the system and the connection relationships between commodities and technologies.

The upper part of Figure 3.1 displays the RES diagram of a small illustrative example corresponding to a simple co-generation system. The system consumes natural gas to produce heating power and electricity. The heating power is provided to the buildings connected to the LMES and the generated electricity is sold to the national grid. In this example, the resource commodity is the natural gas, which is noted as GAS in the figure. The supply resources are the heating power provided to the customers and the generated electricity, which are noted respectively as HEAT2 and ELEC. The hot water stored in the thermal storage tank, noted as HEAT1, is an intermediate commodity. The available technologies are CHP units, boilers and thermal storage. The CHP units consume natural gas to simultaneously produce heat and electricity. The produced heat can be either stored as hot water in the storage tank or transferred directly to customers' end. The electricity generated is sold to the national grid. The hot water produced by boilers is used in the same way as the one produced by CHP components.

However, a RES diagram only displays the technologies that may be used in the system but does not provide the detailed catalog describing the types of device available for each technology. It is thus possible to introduce a more detailed description of the design problem through the use of a diagram representing the superstructure of an LMES. As defined by [START_REF] Yokoyama | Optimal design of gas turbine cogeneration plants in consideration of discreteness of equipment capabilities[END_REF], the superstructure of an LMES is composed of all the units of equipment considered as candidates for selection. The superstructure thus describes the set of available conversion and storage technologies with their corresponding input and output commodities together with a list of the available equipment types for each technology (described e.g. by their size or capacity, their energy conversion performance and their price).

The bottom part of Figure 3.1 displays the superstructure of our illustrative example. Let us assume that, for the design of this co-generation system, two types of CHP units with different efficiency coefficients and capacities are available in the catalog and that we can install at most 5 units for each type in the system. There are also two types of boilers with distinct efficiencies and capacities and we are allowed to install no more than 10 units of each type of boiler. Finally, there is a single thermal storage tank and its storage capacity should not exceed 50 MWh. In the superstructure provided in Figure 3.1, the input parameters associated to the maximum number of each type of devices or to the maximum storage capacity allowed are displayed in grey.

Designing an LMES will thus consist in determining the number of units of each type of devices (among the ones belonging to the superstructure) to be installed in the actual system and in fixing the capacity of the energy storage devices to be installed.

In our simple example, the design decisions would consist in fixing the number of CHP units and boilers of each available type to install (if any) and in determining the size of the thermal storage tank to be built (if any). These decisions are denoted by blue question marks in the superstructure provided in Figure 3.1.

Design considerations

Besides the description of the superstructure, some other aspects need to be taken into account when designing an LMES.

Price and availability of the resource commodities

As mentioned in Subsection 3.2.1, the resource commodities are often bought from an external supplier. The corresponding unit buying price is likely to vary over the lifetime of the LMES. In fact, depending on the case, this price may vary over the day, the week and the season and may also display long-term trends over the years. When designing an LMES, it is important to take into account these variations and thus to follow as precisely as possible the time at which the resource commodities are used. This will namely directly impact our ability at accurately estimating the operation cost of the LMES. Moreover, the daily variations of the price will have a strong impact on the sizing of the thermal storage. In case there are strong variations in the buying price within the day, it may be interesting to install a rather large thermal storage capacity. This will enable us to produce energy in periods where the price of the resource commodities is low, store it for a few hours and use it afterwards to meet the energy demand in periods where the price of the resource commodities is higher. On the contrary, a thermal storage might not be as useful in case the price is relatively stable throughout the day.

Furthermore, we should also take into account the limited availability of these resource commodities. For instance, the electricity, when considered as a resource commodity, is bought from an external supplier through the national electric grid. There is usually a limit on the electric instantaneous power that may be consumed by the system. The value of this limit is negotiated each year with the supplier and gives rise to a contract. The corresponding cost is proportional to the value of the CHP STO Boiler

G A S H E A T 1 H E A T 2 E L E C

Type1: ?/5 contracted maximum electric power. As for renewable energy sources such as solar radiation, as mentioned above, their availability depends on weather conditions and is also strongly limited.

Demand for the supply commodities

The energy demand to be satisfied by an LMES depends on many factors, among which is the ambient temperature, the number of consumers located in the district and their daily and weekly schedules. As a consequence, it displays a strong temporal variability. First, the energy demand strongly varies throughout the day. For instance, the demand for cooling power is usually much higher at daytime than at night. Second, the demand pattern of a weekday often differs from the one observed during the week-end due to the fact that consumers do not have the same routines in working days as in the week-end. Third, the cooling and heating demands are highly dependent on the ambient temperature and thus feature strong seasonal variations: the total cooling (resp. heating) demand of a winter day is thus significantly smaller (resp. higher) than the one of a summer day. These demand variations display a seasonal and periodic pattern and are thus to some extent predictable. Nevertheless, extreme weather conditions which are more difficult to accurately anticipate may happen. These extreme conditions lead to the presence of extreme days, i.e days in which the energy demand is exceptionally high or low as compared to its average value. A major requirement is that the LMES should be able to meet this demand at all time, whatever the hour of the day, the day in the week and the time in the year. In particular, it should be able to satisfy all the demand, even during the extreme days.

Furthermore, an LMES is a long-term investment with an expected lifetime spanning several decades and is often built together with the district it will serve. This means that the forecasted yearly energy demands feature an upward trend over the years as the district is developed and new buildings progressively connect to the LMES. As a consequence, the design of an LMES should not be seen as a one-time decision but rather as a process in which investment decisions are made step by step, following the development of the district and the increase of the demand. This implies that we not only have to determine the final optimal layout of the system but also a multi-phase strategic deployment plan describing how the various conversion and storage devices should be progressively added to the system over the years.

Ambient temperature

The ambient temperature has an impact on the efficiency of some technologies. For example, the COP of electric chillers varies with the ambient temperature and this impact should be considered when estimating the operation cost of the system.

Maintenance cost

Each device installed in an LMES needs to regularly undergo maintenance to stay in a good working condition. The total maintenance cost of a device over its lifetime can be straightforwardly computed once we know when it will be installed in the system. Depending on the project, the maintenance cost of a device is either fixed or varies with the number of years the device functions in the LMES. In both cases, this maintenance cost can be incorporated in the fixed installation cost, resulting either in a fixed or a time-varying total installation and maintenance cost.

Discount rate

Since the construction of an LMES is a long-term investment project lasting several decades, we need to take into account the time value of money. Basically, investing money in the project now is more costly than investing money in the future. This is because an amount of capital, if not invested in the project now, will grow each year thanks to e.g. bank interests.

Let β denote the discount rate, i.e. the percentage by which an amount of money will grow each year if not invested in the project. Investing an amount R of money n years from now is considered equivalent to investing the amount of money R (1+β) n now. Consequently, installing a device in the LMES in the future costs less than installing it now.

Our objective in this work will thus be to minimize the net present cost of the LMES, i.e. the net present value (NPV) of the total cost of designing and operating the system over its whole lifetime.

Problem statement

In Sections 3.2 and 3.3, we described the various elements to be considered when optimally designing an LMES. This optimization problem can be formally stated as follows. Given the available energy resources, the predicted total demand for supply commodities in the neighbourhood over the LMES lifetime and the predefined superstructure, select the technologies to be installed in the system and determine the number and size (i.e. capacity) of the corresponding conversion and storage devices so as to minimize the total construction and operation cost of the system while satisfying the customers' demand at all time.

We now describe in more details the problem input data, the decisions to be made, the main constraints to be complied with and the objective to be reached.

Input data

The problem input data first consist in a description of the superstructure of the system. As explained in Subsection 3.2.3, the superstructure gives a list of all the conversion and storage devices considered as candidates for selection. Each candidate device corresponds to a given technology (with its input and output commodities), to a size (or maximum capacity), to an energy conversion performance (if relevant) and to an installation and maintenance price. We also know the unit cost per kilowatt (KW) of the maximum electric power to be contracted with the electricity supplier.

Second, we are given a set of time series describing the forecasted time variations (usually on a hourly basis) of the buying price and availability of the resource commodities, of the demand for supply commodities and of the ambient temperature over the LMES lifetime.

Finally, the value of the rate to be used to discount the value of future costs is also assumed to be known.

Decisions

Design decisions The design decisions to be made consist in selecting, among the list of technologies and devices provided in the superstructure, the ones to be installed in the system and in determining the value of the maximum electric power to be contracted with the electricity supplier. Note that, regarding the sizing of the installed devices, two different situations occur. For some technologies such as thermal storage, the capacity of the device can be fixed to any continuous value in a predefined range and we only have to determine this value. In contrast, for other technologies such as chillers, CHP units or heat pumps, the capacity of the devices is to be chosen within a discrete list of values corresponding to off-the-shelf models available at the equipment manufacturers. In this case, we should determine the discrete number of devices of each available type to be installed in the system.

Moreover, as mentioned above, as the energy demand features long-term trends over the years, these decisions should not be made once and for all at the beginning of the system lifetime, but rather progressively, year after year, to adjust to the longterm variations of the energy demand. In other words, we should be able to provide, not only a final layout of the system, but also a multi-phase deployment plan covering the system lifetime and describing how the devices should be installed year after year.

Operation decisions Our design decisions will have a strong impact on the daily operational management of the system, in particular on our ability at satisfying the demand for supply commodities. Furthermore, the operation cost, which corresponds mainly to the cost of buying the resource commodities, makes up for a large part of the total cost of the system (typical value is between 20% and 60% of the total cost) and can therefore not be neglected when computing the total cost of an LMES.

The main challenge here is that it is not possible to use average values of the demand for supply commodities and/or average values of the buying price of resource commodities to estimate the impact of our design decisions on the operational management of the system. Namely, both the demands and prices display strong time variations. Moreover, the technical features of the conversion devices, in particular the fact that they have a non-zero minimum working load and a maximum working load, make it impossible to ensure that an LMES able to satisfy the average hourly demand will be able to actually satisfy the time-varying demand at all time. Furthermore, some devices such as chillers have a load-dependent and temperature-dependent COP, which means that the consumption of resource commodities is not proportional to the production of supply commodities. Finally, as the buying price of resource commodities is time varying, it is not enough to know how much resource commodities we need: we also have to determine when these resource commodities will be bought to accurately compute the corresponding cost. Hence, using an average or aggregated representation of the demand to estimate the operation cost may lead to a significant underestimation of its actual value.

As a consequence, in order to accurately estimate the impact of our design decisions on the operational management of the LMES, we should build a detailed operational schedule describing, on a hourly basis, the on/off status and the load allocation of each conversion device, the amount of energy stored in each storage device and the amount of supply commodities bought from an external supplier. Ideally, we should build such a schedule for each hour, each day and each year of the LMES lifetime.

The design and operation scheduling decisions to be made for the illustrative simple cogeneration system are partially represented in Figure 3.2. The upper part of the figure shows the design decisions which have been made for a given year of the system lifetime: note how the blue question marks of the superstructure of Figure 3.1 are replaced by blue numbers. Thus, at the beginning of the considered year, the LMES comprises 2 CHP unites of type 2, 8 boilers of type 1, 2 boilers of type 2 and a thermal storage capacity of 30 MWh. For each hour of the corresponding year, we should build a schedule describing among others the on/off status of each devices. The bottom part of Figure 3

Constraints

The designed LMES should comply with some restrictions. For instance, we have an upper bound on the number of devices of a given type that may be installed. For technologies such as thermal storage for which the capacity can be customized, we also have an upper bound on the capacity of the tank to be installed. We also have a large set of constraints to ensure that the operational schedule built for each hour is relevant in practice. These constraints guarantee among others that the demand for supply commodities is satisfied, that the energy balance in thermal storage devices is respected, that energy conversion is carried out in each device following its actual COP... and this at all time, i.e. at each hour of the LMES lifetime.

Objective

As mentioned above, our objective is to minimize the net present total cost of the LMES, i.e. the net present value (NPV) of the total cost of designing and operating the system over its whole lifetime.

Case studies

We now introduce three case studies which we will use in our numerical experiments to evaluate our models and algorithms. These case studies correspond to LMESs that EDF is currently developing in China.

District cooling systems

A district cooling system (DCS) is a local energy system which distributes cooling capacity in the form of chilled water to a set of buildings located nearby. The customers purchase the delivered chilled water from the DCS operator and use it for air conditioning and dehumidifying or for industrial process cooling.

In a DCS, we have a single supply commodity, the chilled water, whose demand displays strong daily, weekly and seasonal variations. In our numerical experiments, we consider two DCSs located in two different cities. City A is in a tropical region in the south of China so that customers need cooling power during most of the year. City B is in a region with a temperate climate: hence, the cooling demand is mainly concentrated around the summer. In both cases, the expected lifetime of the DCS is 30 years. The population of the district, and as a consequence the demand for chilled water, is anticipated to increase during the first years and to stay stable afterwards.

In the considered DCSs, there is a single resource commodity corresponding the electricity bought from the national grid. The price for electricity is daily periodic but is otherwise stationary and does not change from one day to the next. Its availability is limited by a contract with the electricity supplier that may be renegotiated each year.

In terms of technologies, the superstructure comprises a set of single-mode electric chillers, dual-mode electric chillers and a thermal storage. The capacity of the thermal storage tank may take any continuous value within a predefined range. As for the chillers, we have a discrete set of available models that may be bought from equipment manufacturers and we should determine the number of chillers of each type to be installed in the DCS. Figure 3.3 displays the RES diagram of the DCSs considered in our two case studies. The resource commodity is denoted by ELEC, the intermediate commodity by ICE and the supply commodity by COLD. Single-mode chillers consume the ELEC commodity to produce COLD. As for the dual-mode chillers consume, they also consume the ELEC commodity but produce both ICE to be stored in the storage tank and COLD to be directly delivered to the customers' buildings.

Trigeneration system

The third case study corresponds to a trigeneration system located in a city where both heating and cooling power are demanded. The system serves various types of customers in the surrounding area, such as residential buildings, commercial centers, office buildings and hotels. Its expected lifetime is 30 years. The served district is currently under development: the population is thus expected to increase for the first four years and to stay stable afterwards.

The system should satisfy both a cooling demand and a heating demand. These demands have a seasonal feature. Cooling power is needed during 6 months whereas heating power is required during 4 months. The cooling supply period and heating supply period have no overlap. We assume that the energy demand grows with the same trend as the population in the district. The ambient temperature of the system is considered to be yearly periodic.

The available resource commodities are natural gas and electricity. The unit buying price natural gas is constant in time. The price for electricity is daily periodic but is otherwise stationary and does not change from one day to the next.

The technologies available in the superstructure include CCHP units, boilers, air source heat pumps, single-mode electric chillers, dual-mode electric chillers, hot water thermal storage and ice thermal storage. Recall that for thermal storage, the capacity can be set to any value in a predefined continuous range. For the other technologies, the installed devices should be selected among the ones available in the manufacturer's catalog so that there is only a finite and discrete set of potential capacities.

Figure 3.4 displays the RES diagram of the trigeneration system. On this figure, ELEC1 and GAS denote the resource commodities, ELEC2, HEAT1 and ICE the intermediate commodities and HEAT2 and COLD the supply commodities. The CCHP unit converts the resource commodity GAS into electricity ELEC2 (to be consumed within the LMES), heating power HEAT1 (to be stored in the thermal storage STO_HEAT), HEAT2 and COLD (to be distributed directly to the customers' buildings). As for the boiler, it also converts GAS into HEAT1 or HEAT2 heating power. We then have three technologies powered by electricity (either bought from the supplier ELEC1 or generated locally by the CCHP unit ELEC2). The heat pump produces HEAT1 or HEAT2 heating power and COLD, i.e. cooling power to be directly distributed to the customers in the form of chilled water. Single-mode chillers only produce COLD while dual-mode chillers produce both COLD and ICE, i.e. cooling power to be stored in the thermal storage STO_COLD before being distributed.

Conclusion

This chapter provided a detailed description of the optimization problem under study in this PhD work. We first introduced the main components of an LMES by focusing on the ones which have a major impact on the total cost of the system. These components are roughly classified into two classes: commodities and technologies. This classification allowed us to represent the layout of an LMES using either an RES diagram or a more detailed superstructure diagram. Second, we presented several additional key aspects of the problem which also have a strong influence on the LMES final layout. The overall optimal design problem was then formally stated as a classic combinatorial optimization problem. We underlined the necessary input data, the decision variables, the main constraints to be respected and the objective function. Finally, we briefly introduced the three case studies that we will use in our numerical experiments. They correspond to two types of energy systems: district cooling systems and trigeneration systems.

Chapter 4

Problem modeling and mathematical formulation

Introduction

The combinatorial optimization problem described in Chapter 3 leads to the formulation of a huge-size mixed-integer non-linear program which is computationally intractable. This huge size mainly comes from the need to simultaneously make strategic design decisions over a long-term planning horizon and evaluate the impact of these strategic decisions on the operation cost of the system through the building of detailed schedules using short (typically hourly) time-steps. Moreover, the presence of energy conversion devices whose efficiency depends on the load (e.g. electric chillers) brings some non-linearities in the mathematical programming formulation.

In this chapter, we first present in Section 4.2 the modeling choices and assumptions we use to reduce to some extent the size of the obtained mathematical program and to eliminate the aforementioned non-linearities. We then provide in Section 4.3 the resulting MILP formulation. Finally, Section 4.4 presents the computational results obtained when directly solving this MILP formulation with a mathematical programming solver on instances based on the three case studies described in Section 3.5.

Problem modeling

Investment phases

The lifetime of an LMES usually spans several decades: let us denote by Y the corresponding number of years. As explained earlier in Section 3.4, the energy demand to be satisfied by the LMES is likely to display a long-term upward trend, especially when the LMES is built in a new district currently under development. This means that it may not be cost-efficient to build the whole system in one step. We should rather only install part of the energy conversion and storage devices at the beginning of the system lifetime and progressively add, year after year, the other devices as they become necessary to adjust to the long-term variations of the demand.

In what follows, we propose to divide the system lifetime into a set of Φ investment periods termed phases. A phase typically spans one or several years. Let y ϕ (resp. y ϕ ) be the index of the first (resp. the last) year belonging to phase ϕ.

We make two modeling assumptions here. First, design decisions, such as the installation of new devices or the building of additional storage capacity, can only be made at the beginning of a phase. Second, the hourly evolution over the year of the energy demands, of the price and availability of resource commodities are described by the same time series for all the years belonging to the same phase. Note that these two modeling assumptions are based on the ones currently used by EDF practitioners when studying the design of an LMES.

Representative days

Ideally, the operation cost during a given investment phase should be estimated by building a detailed operation schedule describing the on/off status and the load allocation of each device together with the stored energy for each hour, each day and each year of this phase. This clearly leads to a huge intractable optimization problem.

A first simplification can be made here by exploiting the second assumption on investment phases mentioned in Subsection 4.2.1. Basically, this one amounts to assuming that all years belonging to the same phase are strictly identical in terms of energy demand, commodity prices and availability. Thus, a single operation schedule covering 365 days may be built for all these years and the total operation cost over the investment phase can be computed as the sum, over all years belonging to this phase, of the actualized value of the cost of this schedule.

However, building an operation schedule covering the 365 days of a full year for each investment phase is still out of reach in terms of computational effort. This is why, in order to further decrease the size of the problem, we propose to build an operation schedule only for a subset of disjoint typical and extreme days. These days should be carefully selected in order to represent, as best as possible, the various conditions under which the LMES will be operated and to obtain an estimation, as accurate as possible, of the operation cost.

We propose to select a subset of typical days based on an analysis of the time series representing the hourly evolution of the input parameters over a full year. More precisely, we seek to identify the 'best' typical days by solving a clustering problem which corresponds to an extension of the k-medoid problem investigated by [START_REF] Zatti | k-MILP: A novel clustering approach to select typical and extreme days for multi-energy systems design optimization[END_REF]. This problem consists in partitioning the set of days belonging to the available time series into a predefined number of groups or clusters and in selecting one member (i.e. one day) in each cluster to represent all the days belonging to the cluster. The objective is to minimize the sum of the Euclidean distance between each day in the time series and its representative.

Moreover, we propose to choose 4 extreme days for each supply commodity and for each resource commodity with a limited availability. For each supply commodity, we first select 2 days, the day with the highest total daily demand and the one with the highest hourly demand, in order to ensure that the capacity of the installed devices will be large enough to satisfy the demand even if it is extremely high. Moreover, in case the demand is very low but not equal to zero, we may encounter difficulties to satisfy it due to the constraints on the minimum production power of some devices. In order to avoid it, we also select as extreme days the day with the lowest non-zero total demand and the one with the lowest non-zero hourly demand. As for the resource commodities with a limited availability (e.g. solar radiation), we also select 4 extreme days: the day with the highest hourly availability, the one with the highest total daily availability, the one with the smallest hourly availability and the one with the smallest total daily availability. Note that, the extreme days are assumed to happen once in the year so that the corresponding cost has a small impact on the total annual operation cost. However, by imposing that the system is able to satisfy the energy demands even under extreme conditions these extreme days have a huge impact on the design of the system. Moreover, for the time being, we do not select extreme days with respect to the price of resource commodities because in our case studies, these prices are either constant or follow the exact same pattern every day.

Figure 4.1 shows the hourly evolution of the cooling demand for 6 representative days which may be selected with respect to the supply commodity COLD to estimate the operation cost during a given investment phase for the district cooling system located in City B (see Subsection 3.5.1). Note that the hourly demand is normalized, i.e., is expressed as a proportion of the maximum hourly cooling demand. We choose 2 typical days and 4 extreme days to represent the cooling demand variation of the year. Each representative day is allocated a weight, which corresponds to the number of days of the year it stands for in the clustering. Thus, Day 136 (resp. Day 263) is a typical day which acts as the representative of a cluster comprising 66 (resp. 95) days and will be assigned a weight equal to 66 (resp. 95). Each extreme day is assumed to occur only once a year and is thus assigned a weight equal to 1. Note that the sum of these weights is not equal to 365 because in the considered DCS located in a region with a temperate climate, there is no cooling demand in more than half of the year.

Piece-wise linear approximation of non-linear performance curves

As explained in Chapter 3, for some technologies, the conversion efficiency coefficient is not constant but rather depends on the ambient temperature and/or on the load. In this case, the manufacturer of each conversion device provides us with a set of non-linear performance curves giving the energy consumption as a function of its load rate under a discrete set of values of the ambient temperature. A typical example are the electric chillers. The conversion efficiency of a single-mode chiller is described by a set of performance curves: each curve depicts the electric power consumed by the chiller as a function of the produced cooling power under a certain ambient temperature. Similarly, the conversion efficiency of a dual-mode chiller is described by two sets of performance curves: one set for its cooling mode and one set for its ice-production mode.

The non-linearity of these performance curves poses a major difficulty for the efficient resolution of the optimization problem as this one would be formulated as a mixed-integer non-linear program.

To overcome this difficulty, we propose, for each type of conversion device featuring a part-load efficiency, to approximate each of its non-linear performance curves through a piece-wise linear function. This function is built using the following procedure. We first choose a number of breakpoints to be used in the piece-wise linear approximation. We then determine the coordinates of these breakpoints so as to minimize the distance between the approximate performance curve and the actual one. This is done by heuristically solving a small non-linear unconstrained optimization problem through the Powell's conjugate direction method [START_REF] Powell | An efficient method for finding the minimum of a function of several variables without calculating derivatives[END_REF]. In Figure 4.2, we show as an example the piece-wise linear approximation of a the performance curve of a single-mode electric chiller (SMEC). Four breakpoints are used to approximate the non-linear curve so that the piece-wise linear approximation comprises three continuous segments. We then consider historical data about the average ambient temperature at each time-step of each representative day and use this information to identify which approximate performance curve (among the ones built in the previous step) should be used to compute the energy consumption of each type of device during this timestep.

Leveraging the convexity of the performance curves

When applicable, we propose to exploit the convexity of the performance curves of the conversion devices to further reduce the size of the optimization problem. When the performance curves of a given type of energy conversion devices are convex, we have the following property. Lemma 4.2.1. For a set of identical active conversion devices producing the same commodity, the optimal load allocation consists in equally distributing the total output power between the active devices.

Proof. The proof is done by contradiction. Suppose we have two identical active devices, producing a total output power of P with device 1 producing P 1 and device 2 producing P 2 > P 1 . Let π : P → Q = π(P) be the convex performance curve of these two devices. The total amount of input commodity consumed by the two devices producing P is Q = π(P 1 ) + π(P 2 ).

We show that this load allocation is not optimal, i.e. that it is possible to reduce the total amount of consumed commodity. Namely, let δP be a small variation in the output. By decreasing the output of device 2 by δP, we can obtain a decrease in the input commodity consumed by this device of π ′ (P 2 )δP, where π ′ is the derivative function of π. In order to still be able to provide a total output of P, we increase the output of device 1 by δP, which leads to an increase in its input commodity consumption of π ′ (P 1 )δP.

By convexity of function π, π ′ (P 1 ) < π ′ (P 2 ). Hence the total amount of input commodity consumed with the load allocation (P 1 + δP, P 2 -δP) is smaller than the one consumed with the load allocation (P 1 , P 2 ). The result follows.

Let us now focus on the case in which the performance curve π is convex and piece-wise linear. When several identical active devices are simultaneously producing the same output commodity, the relation providing the total amount of consumed input commodity as a function of the total amount of output power can be plotted as an aggregate performance curve. We have the following property: Lemma 4.2.2. Let π be the piece-wise linear and convex performance curve of a given type of active device producing a given commodity. π involves B breakpoints. Let (a b , o b ) be the abscissa and ordinate of breakpoint b.

The aggregate performance curve Π S of S identical active devices of this type simultaneously producing the same commodity is also piece-wise linear and convex. It involves B breakpoints whose coordinates are given by (Sa b , So b ).

Proof. Let us consider the case where the total output of the S active devices is P ∈ [SP min ; SP max ] where [P min , P max ] is the production range of a single device.

By Lemma 4.2.1, the optimal load allocation consists in requiring that each active device σ = 1...S produces the same output P By generalizing this result to all possible values of the total output P, we have that Π S is a piece-wise linear function involving B breakpoints of coordinates (Sa b , So b ). Moreover the slope of Π S on its b th segment is s b . As π is convex, we have

s b ≤ s b+1 , b = 1...B. Π S is thus convex.
We first note that, when the performance curves are convex, the optimal load allocation between a set of identical active devices producing the same commodity consists in equally distributing the total output power between the active devices: see Lemma 4.2.1. For instance, a set of identical single-mode electric chillers simultaneously producing cooling power should be operated at the same working load to minimize the electricity consumption.

Moreover, let us now focus on the case where the performance curve of these devices is piece-wise linear and convex. In this case, when several identical active devices are simultaneously producing the same output commodity, the relation providing the total amount of consumed input commodity as a function of the total amount of output power can be plotted as an aggregate performance curve taking the form of a piece-wise linear and convex function: see Lemma 4.2.2. More precisely, assume that the performance curve of an individual active device producing a given commodity involves B breakpoints with (a b , o b ) the abscissa and ordinate of breakpoint b. The aggregate performance curve of S identical active devices of this type simultaneously producing the same commodity is also piece-wise linear and convex and involves B breakpoints whose coordinates are given by (Sa b , So b ). Figure 4.3 illustrates this for a given type of single-mode electric chillers. The blue curve describes the performance of a single active chiller of this type as a piece-wise linear function. The orange (resp. green) curve describes the aggregate performance of two (resp. three) chillers of this type simultaneously producing cooling power at the same load rate. Thus, each point of these aggregate performance curves gives (in ordinate) the minimum electricity power needed to produce the demanded cooling power (in abscissa) for the corresponding number of identical active chillers. Note that in case the performance curves are linear, i.e. in case the energy conversion coefficient is constant, any load allocation of the total output power of a set of identical active devices leads to the same total consumption of the input commodity. It is thus also possible to build aggregate linear performance curves in this case.

The use of these (convex or linear) aggregate performance curves enables us to simplify to some extent the building of the operation schedules for each representative day. Namely, it is possible to build these schedules by using aggregate scheduling variables describing the number of on/off devices of each type, their aggregate production of output commodity and their aggregate consumption of input commodity (rather than detailed scheduling variables describing the on/off status, output and input power of each individual conversion device). This leads to a significant decrease in the number of variables and constraints introduced in the formulation. It also avoids symmetry problems coming from the fact that there exist many alternative optimal solutions that differ from one another only by the way each detailed operation schedule is assigned to one of the identical conversion devices.

Moreover, thanks to the fact that we will seek to minimize the operation cost, we do not have to introduce extra binary variables to handle these piece-wise linear and convex functions. We only have to add a set of linear constraints to ensure that the point representing the aggregate amounts of output and input commodity of a set of S identical active devices will be in the epigraph of each corresponding aggregate curve. In any optimal solution, this point will be on the aggregate performance curve. This leads to the formulation of the union of a sequence of disjoint convex sets depending on three aspects: the aggregate cooling power, the aggregate electricity consumption power and the number of active identical chillers. Such formulation problem is usually called a disjunctive programming problem and we recommend readers to refer to [START_REF] Pablo | Mixed integer linear programming formulation techniques[END_REF] for more details.

Mixed-integer linear programming formulation

Thanks to the modeling choices and assumptions introduced in the previous section, the problem of optimally designing an LMES can be formulated as an MILP. We provide in this section a detailed description of the corresponding formulation.

Indices and input parameters Time representation and management

We use a three-level time discretization in our problem modeling. We first divide the system lifetime involving Y years into a set of Φ investment periods termed phases. Let y ϕ (resp. y ϕ ) be the index of the first (resp. the last) year belonging to phase ϕ. The index d represents the time unit used to build schedules at the operational level. In practice, the length of this time unit is chosen from a half day to one week. In this thesis, we use a time unit of one day to build operation schedules. For each considered phase ϕ ∈ {1, ..., Φ}, a set D ϕ of typical and extreme days selected as explained in Subsection 4.2.2 are used to represent the various daily patterns of the energy demand, commodity price and availability. Each representative day d of phase ϕ, denoted by (ϕ, d), is assigned a weight w ϕ,d . If d is a typical day identified through the clustering approach discussed in Subsection 4.2.2, w ϕ,d corresponds to the number of days it represents, i.e. to the number of days in the original data assigned to the same cluster as day (ϕ, d). If d is an extreme day, w ϕ,d is set to 1. Finally, in order to depict the intra-day variation of the input parameters, each representative day is further divided into H time-steps indexed from 0 to H -1, each of which has a duration of ∆t hours. For the sake of readability, in the remainder of this thesis, we use the bold letter t to represent the time period (ϕ, d, h), which corresponds to phase ϕ, selected day d and time step h.

The net present value of the future cost is calculated using an annual discount rate denoted by γ. We use two distinct actualization rates in our problem formulation. The first one, called the design actualization rate, is used to compute the net present value of the system design cost at the beginning of each phase. It is equal to α ϕ = 1 (1+γ) y ϕ for phase ϕ. The second one, called the operation actualization rate, is used to calculate the net present value of the operation cost during phase ϕ. It is defined by β ϕ = ∑ y ϕ y=y ϕ 1 (1+γ) y for phase ϕ. Recall that all years belonging to a given phase ϕ are assumed to be described through the same input data so that the total annual operation cost is constant for all years belonging to ϕ. The net present value of the operation cost during ϕ can thus be obtained by multiplying the total annual operation cost during phase ϕ by β ϕ .

Commodities

Let N = {GAS, ELEC, HEAT, COLD, ICE, ...} denote the various forms of energy in the considered LMES. As can be seen from the examples and case studies described in Chapter 3, in an LMES, a given form of energy may correspond to different commodities ( e.g. electricity bought from the national grid or generated by a gas engine). In order to distinguish between them, we denote each commodity by a pair (n, s): index n ∈ N represents the form of energy and index s ∈ {1, • • • , S n } is used to number commodities with the same type n. The set of all commodities is thus C = {(n, s) : n ∈ N , s ∈ {1, • • • , S n }}. C may be partitioned into the set of resource commodities C R , the set of intermediate commodities C I and the set of supply commodity C S .

For example, the set of commodities involved in the trigeneration system described by Figure 3.4 is given by: C = {(ELEC, 1), (ELEC, 2), (GAS, 1), (HEAT, 1), (HEAT, 2), (ICE, 1), (COLD, 1)} with C R = {(ELEC, 1), (GAS, 1)}, C I = {(ELEC, 2), (HEAT, 1), (ICE, 1)} and C S = {(HEAT, 2), (COLD, 1)}.

We use EP t,c to denote the energy price of a commodity c at timestep t. This parameter corresponds to the price for purchasing one unit of energy if c ∈ C R is a resource commodity. Moreover, let C S,Sell ⊂ C S be the subset of supply commodities which may be sold to an outsider buyer in addition (or instead of) being sold to the buildings connected to the LMES (e.g. electricity generated by the LMES sold to the national grid). For c ∈ C S,Sell , EP t,c denotes the price for selling one unit of energy to the outsider buyer.

The demand of a supply commodity c at time-step t, i.e. the amount of this commodity to be produced by the LMES, is denoted by Dem t,c . In a typical cooling system, the difference between the temperature of the circulating chilled water and the ambient temperature is less than 20 • C. Therefore, we may assume that the thermal loss during distribution in a cooling system is negligible and directly use the actual demand of the customers for cooling power to set the value of Dem t,c . In contrast, in a heating system, the difference between the temperature of the circulating hot water and the ambient temperature is above 60 • C. We should thus consider heat losses during distribution. The heat loss rate, denoted by ϵ t,HEAT , is the ratio of heat dissipation when the hot water is transferred from the LMES to customers at time-step t. To take this into account, Dem t,c should be set to the value of the actual demand of the customers for heating power divided by (1ϵ t,HEAT ).

Finally, the availability of some resource commodities is limited. Let C R,Ph be the set of resource commodities for which this limited availability comes from physical reasons (e.g. renewable energy resources such as solar radiation whose availability depends on weather conditions). In this case, the amount of commodity c available at time t is an input parameter denoted by Conso max t,c . Moreover, the availability of some resource commodities (e.g. electricity bought from the national grid) may also be limited due to a contract signed with the outside supplier. Let C R,Co be the set of resource commodities for which this limited availability comes from a contract negotiated with an outside provider. Let MaxC c be the highest allowed value for this contracted maximum instantaneous power. In theory, the contracted maximum power may take any value between 0 and MaxC c . However, in practice, this one is most often fixed at a value corresponding to an integer multiple of a predefined discretization unit denoted by C

Technologies

Let P = {CCHP, BOILER, SMEC, ..., STO_HEAT, ...} denote the set of available energy conversion and storage technologies. Each technology is described among other by N in p , the set of energy forms it consumes, and N out p , the set of energy forms it produces. In what follows, for the sake of simplicity, we will assume that |N in p | = 1 and denote by n in p the single form of energy consumed by technology p. By using the high-level description of the LMES structure provided by the RES diagram, we may define, for each technology p ∈ P, the set C in n,p (resp. C out n,p ) of all consumed (resp. produced) commodities belonging to the energy form n ∈ N . The other way round, we may define, for each commodity c ∈ C, the set P out f low c (resp. P in f low c

) of technologies that consume (resp. produce) commodity c. As explained in Chapter 3, the set of technologies P can be partitioned according to the nature of the design decisions to be made. Thus, for the set P DIS ⊂ P of discrete technologies, the design decisions consist in selecting from a discrete list of equipment types available in the manufacturers' catalogs the number of devices of each type to install. In this case, let l ∈ {1, • • • , L p } denote the various types of device available for technology p. The pair m = (p, l) thus represents a given type of energy conversion device belonging to technology p. The maximum number of devices of type m that may be installed is denoted by MaxNum m .

In contrast, for the set P CONT ⊂ P of continuous technologies, the design decisions consist in selecting the capacity of the single piece of equipment to be installed within 0 and a maximum size. In this second case, for the sake of readability, we will assume L p = 1 and denote by (p, 1) the single available type of equipment, the capacity of which should be determined during the design of the system. The maximum size of the installed piece of equipment is denoted by MaxSize m . Although in theory, the size of the installed device may take any value between 0 and MaxSize m , in practice, this one is most often fixed at a value corresponding to an integer multiple of a predefined discretization unit denoted by Size step m . The maximum number of discrete units of capacity that may be installed for the device of type m is thus

MaxNum m = ⌊ MaxSize m Size Step m ⌋. Let M = {(p, l)|p ∈ P, l ∈ {1, • • • , L p
}} be the set of all types of devices and M p 1 = {(p, l) ∈ M|p = p 1 } be the subset of all devices belonging to technology p 1 ∈ P.

Each type of device m ∈ M is associated with an installation cost FC m and an accumulated maintenance cost MC ϕ,m . The installation cost includes the purchase and construction costs, to be paid only once at the beginning of the phase in which the device is installed. Note that FC m corresponds to the cost of installing one device m = (p, l) if m belongs to a discrete technology p ∈ P DIS and to the cost of installing one discrete capacity unit of size Size step m if m belongs to a continuous technology p ∈ P CONT . Once this device (or discrete capacity unit) is installed, we may compute its actualized total maintenance cost over the interval between its installation time (beginning of phase ϕ) and the end of the system' lifetime. The corresponding cost, whose value depends on the installation phase, is defined as the accumulated maintenance cost MC ϕ,m . Therefore, when installing a device of type m in phase ϕ, the net present value (NPV) computed at phase ϕ of the total cost is FC m + MC ϕ,m and the NPV computed at the beginning of the project can be obtained by multiplying FC m + MC ϕ,m by the design discount rate α ϕ .

In what follows, we will also distinguish when needed between the set of conversion technologies P CONV ⊂ P and the set of storage technologies P STO ⊂ P. We denote by M CONV the set of types of energy conversion devices and M STO the set of types of energy storage devices.

Regarding the energy conversion devices, we have a set of additional input parameters to describe the working range and energy conversion efficiency.

We thus define, for each type of conversion device m = (p, l) with p ∈ P CONV and each form of energy n ∈ N out p produced by technology p, P min n,m (resp. P max n,m ) the minimum (resp. maximum) amount of energy of type n that an active device of type m can produce within one scheduling time-step.

As for the conversion efficiency, two main situations may occur:

• The efficiency does not depend on the load. In this case, we use the input data describing the forecasted evolution of the ambient temperature over the year to determine the value of the energy conversion efficiency to be applied to each scheduling time-step. We thus introduce a coefficient η t,n,m giving the amount of energy of form n ∈ C out n,p produced per unit of consumed energy of type n in p by an active device of type m at scheduling time-step t. Note that, when the energy conversion efficiency depends neither on the load nor on the ambient temperature, the value of coefficient η t,n,m is the same for all time-steps t.

• The efficiency depends on the load. In this case, we use the procedure described in Subsection 4.2.3 to build a piece-wise linear approximation of the performance curves provided by the manufacturer and to determine which approximation should be used at each time-step t. Let B t,n,m be the number of breakpoints used in the approximation of the performance curve providing the amount of energy of form n ∈ N out p produced per unit of consumed energy of type n in p by an active device of type m at scheduling time-step t. Each of these breakpoints has coordinates given by (a t,n,m,b , o t,n,m,b ).

Design variables and constraints

In order to formulate the problem as an MILP, we first introduce the decision variables relative to the system layout at each phase ϕ. Note that, to avoid ambiguity in the notation, the names of all decision variables are printed in bold letters.

For each discrete technology p ∈ P DIS , we introduce SD ϕ,m , the integer decision variable representing the number of devices of type m = (p, l) ∈ M p added to the system at the beginning of phase ϕ.

Similarly, for each continuous technologies p ∈ P CONT , we introduce SD ϕ,m , the integer decision variable indicating the number of discrete capacity units of type m = (p, l) ∈ M p added the system at the beginning of phase ϕ.

This leads to the following design constraints, stating that, for each m ∈ M, the total number of devices /discrete capacity units present in the system during the last investment phase Φ should be less that the maximum allowed number.

Φ ∑ φ=1 SD φ,m ≤ MaxNum m ∀m ∈ M (4.1)
Furthermore, we have a second set of design decisions. These decisions are relative to the resource commodities c ∈ C R,Co bought from an outside supplier under a contract specifying a maximum value for the instantaneous power of this commodity that may be bought. This maximum value may be changed from one phase to the next. We thus introduce Cdis ϕ,c the integer decision variable indicating the number of discrete units of power contracted at phase ϕ with the supplier of commodity c.

This leads to the following design constraints, ensuring that, for each commodity c ∈ C R,Co and each phase ϕ, the contract signed with the provider of c stays below the maximum allowed value.

Cdis ϕ,c ≤ MaxUnits c ∀c ∈ C R,Co , ∀ϕ (4.2)
In what follows, we detail the variables and constraints introduced in the formulation in order to build a feasible operation schedule for each day d ∈ D. We first describe the variables and constraints relative to the commodities, then the ones relative to the technologies and finally the constraints linking commodities and technologies.

Operation variables and constraints for commodities

Variables for the commodities

In the operational schedule relative to time-step t = (ϕ, d, h), the total amount of commodity c produced (resp. consumed) by all the devices corresponding to technology p ∈ P is indicated respectively by the continuous variable Inflow t,c,p (resp. Outflow t,c,p ). Note that, thanks to this notation, each arrow in a RES diagram is associated with a flow variable: if the arrow points towards the vertical line representing the commodity, it will be associated to an Inflow variable and if it starts from this vertical line, it will be associated to an Outflow variable. Therefore, at time-step t, the production of commodity c 1 by technologies p 1 and p 2 is denoted as Inflow t,c 1 ,p 1 and Inflow t,c 1 ,p 2 , and the consumption of commodity c 1 by p 3 is denoted as Outflow t,c 1 ,p 3 .

Moreover, we introduce two additional set of decision variables for each timestep: Conso t,c , the total amount of resource commodity c ∈ C R consumed by the system within t, and Sold t,c , the amount of supply commodity c ∈ C S,Sell sold within t.

Table 4.1 recalls the sets of indices, data and variables relative to commodities.

We now describe the constraints relative to commodities to be added to the MILP formulation. These constraints are grouped according to the type (resource, intermediate or supply) of the considered commodity.

Constraints for the resource commodities

The total amount of resource commodity c ∈ C R consumed by the system during time-step t is equal to the sum of the amount consumed by all the technologies in P out f low c .

Conso t,c = ∑ p∈P out f low c Outflow t,c,p ∀c ∈ C R , ∀t (4.3) Set of indices N Set of energy forms C Set of commodities C R Set of resource commodities C R,Ph
Set of resource commodities with a physically limited availability C R,Co

Set of resource commodities with a contractually limited availability

C I Set of intermediate commodities C S Set of supply commodities C S,Sell
Set of supply commodities sold to an outside buyer P Total amount of resource commodity c ∈ C R consumed by the system within t (Continuous)

Sold t,c

Amount of supply commodity c ∈ C S,Sell sold within t (Continuous) For each resource commodity c ∈ C R,Ph , the total amount consumed by the system in each time-step t should stay below the maximum available amount determined by some physical limitations.

Conso t,c ≤ Conso max t,c ∀c ∈ C R,Ph , ∀t (4.4)
For each resource commodity c ∈ C R,Co , the total amount consumed by the system in each time-step t should stay below the maximum available amount determined by the contract signed with the outside supplier. This limit is equal to the maximum instantaneous power contracted for phase ϕ, C step c Cdis ϕ,c , multiplied by the length of the time-step ∆t.

Conso t,c ≤ C step c Cdis ϕ,c ∆t ∀c ∈ C R,Co , ∀t (4.5)

Constraints for the supply commodities

The LMES should be able to satisfy the customers' demand at all time. Thus, for each supply commodity c ∈ C S , the total amount of energy produced by the technologies in P in f low c in time-step t should be equal to the sum of the forecast energy demand coming from the buildings connected to the LMES and of the amount of energy sold to an outside buyer.

The amount of energy sold to an outside buyer is equal to 0 in time-step t for all supply commodities c ∈ C S not belonging to C S,Sell .

Sold t,c =

Operation variables and constraints for the conversion technologies

We now focus on modeling the energy conversion technologies P CONV at the operation level. Recall that M CONV denotes the set of all types of conversion devices: M CONV = {(p, l) : p ∈ P CONV , l = 1, . . . , L p }.

For the sake of simplicity, we will use two assumptions in what follows. First, we suppose that P CONV ⊂ P DIS , i.e. that all conversion technologies are of discrete nature. This means in particular that we assume that there is a set of L p types of devices available for each technology p and that that the total number of devices of type m = (p, l) available at the beginning of phase ϕ can be computed using the discrete design variables as ∑ ϕ φ=1 SD φ,m . Second, we will assume that the performance curves of all available types of conversion devices are either linear (which corresponds to a load-independent conversion efficiency coefficient) or piece-wise linear and convex. As explained in Subsection 4.2.4, this allows us to use aggregate scheduling variables indicating, for a given type of device and a time-step, the number of on/off devices and the aggregate production and consumption of energy and to avoid building a detailed schedule for each individual device.

Variables for the conversion technologies

We first introduce, for each type of energy conversion device m ∈ M CONV , two sets of continuous variables: P in t,n,m (resp. P out t,n,m ) gives the aggregate amount of energy of form n ∈ N in p (resp. n ∈ N out p ) consumed (resp. produced) by the devices of type m in t.

Second, the integer decision variable S t,m represents the number of active devices of type m in t. These variables are needed among others for all devices belonging to technologies for which the minimum output power of an operating device is strictly positive.

Finally, recall that for some technologies, the devices may produce energy under various operating modes: see e.g. CCHP units which may produce either cold or hot water (in addition to electricity) or dual-mode electric chillers (DMEC) which may produce either cold water or ice. These operating modes are mutually exclusive, i.e., in a time-step t, a given device is either off or producing using a single operating mode. Each operating mode will be identified by the distinguishing form of energy it produces. We thus introduce MODE t,n,m the integer variable representing the number of devices of type m operating in mode n ∈ N out p in t. Note that these variables are introduced in the formulation only for the types of devices for which several operating modes exist. Table 4.2 recalls all the notation (set of indices, input parameters and decision variables) related to commodities.

In what follows, for the sake of simplicity, we do not provide operation constraints for all potential technologies but only for the subset of technologies mentioned in the case studies reported in Chapter 3. Note that these constraints are basically of three main types. They ensure that the schedule is feasible with respect to the limited capacity defined by design decisions, that each device is operated within its working range and that the energy consumed by an active device complies with its energy conversion efficiency.

Constraints for Combined Cooling, Heating and Power (CCHP) units

At time-step t, the number of active CCHP units should be less than the total number of CCHP units present in the system, which is computed using the design variables SD ϕ,m .

S t,m ≤ ϕ ∑ φ=1 SD φ,m ∀m ∈ M CCHP , ∀t (4.9)
A CCHP unit converts natural gas either into electricity and heating power or into electricity and cooling power. It thus has two operating modes that may be referred to as heating mode (n = HEAT) and cooling mode (n = COLD). At time-step t, a CCHP unit operates using a single mode. The total number of active CCHP units in t is thus equal to the sum, over the two operating modes, of the number of CCHP units operating in each mode.

A CCHP unit consumes a single type of energy form (n in CCHP = GAS) and its conversion efficiency is load-independent. For each type of devices m ∈ M CCHP , this efficiency is thus characterized by three coefficients: an electricity efficiency coefficient η t,ELEC,m , an heating efficiency coefficient η t,HEAT,m and a cooling efficiency coefficient η t,COLD,m . The relation between the amount of consumed gas and the amount of produced electricity at time-step t is thus given by the following relation. Note that in a CCHP unit, all the waste heat coming from the exhausted gas produced by the internal combustion engine does not have to be transformed into heating or cooling energy. This is why Constraints (4.14) are expressed as inequalities rather than equalities.

P in

Constraints for air source heat pumps (ASHP)

The number of ASHPs active at time-step t should be less than the total number of ASHPs currently installed in the system.

S t,m ≤ ϕ ∑ φ=1 SD φ,m ∀m ∈ M ASHP , ∀t (4.15)
An ASHP converts electricity into either heating power or cooling power. It thus has two operating modes which may be referred to as heating mode n = HEAT or cooling mode n = COLD. The total number of active ASHPs in t is thus equal to the sum, over the two operating modes, of the number of ASHPs operating in each mode.

Mode t,HEAT,m

+ Mode t,COLD,m = S t,m ∀m ∈ M ASHP , ∀t (4.16)
Heat pumps should be operated within their working range. Recall that a heat pump does not have a minimum load rate, i.e. P min n,m = 0 for n ∈ {HEAT, COLD} and m ∈ M ASHP , and can thus be operated at any value between 0 and its maximal capacity. The aggregate amount of energy of form n ∈ {HEAT, COLD} produced by the ASHPs operating in mode n ∈ {HEAT, COLD} in t should thus comply with the following inequality.

P out t,n,m ≤ P max n,m Mode t,n,m ∀m ∈ M ASHP , ∀t, ∀n ∈ {HEAT, COLD} (4.17)
The conversion efficiency of an ASHP does not depend on its load. The relation between the amount of consumed electricity and the amount of produced cooling/heating energy is thus given by the following equality.

P

Constraints for boilers

The number of boilers of type m active at time-step t should be less than the total number of boilers currently installed in the system:

S t,m ≤ ϕ ∑ φ=1 SD φ,m ∀m ∈ M BOILER , ∀t (4.19)
The heating power produced by a boiler can take any value between zero and its maximum capacity, i.e. P min HEAT,m = 0 for any m ∈ M BOILER . At any time-step t, the aggregate production of a set of active boilers of type m should thus lie between zero and the capacity of one boiler multiplied by the number of active boilers.

P out

Constraints for single-mode electric chillers (SMEC)

The number of SMECs of type m active at time-step t should be less than the total number of SMECs currently installed in the system.

S t,m ≤ ϕ ∑ φ=1 SD φ,m ∀m ∈ M SMEC , ∀t (4.22)
A SMEC of type m should be operated within its working range defined by [P min COLD,m ; P max COLD,m ], with P min COLD,m > 0. The aggregate amount of cooling power (n = COLD) produced by the set of SMECs of type m active at t should thus lie within the aggregate working range defined as follows.

P min

COLD,m S t,m ≤P out t,COLD,m ≤ P max COLD,m S t,m ∀m ∈ M SMEC , ∀t (4.23)

The conversion efficiency coefficient of a SMEC depends on its load and on the ambient temperature. We will therefore use the approach described in Subsection 4.2.3 to handle the set of performance curves provided by the manufacturer of each considered type of SMEC. Thus, the conversion efficiency of a SMEC of type

m

Constraints for dual-mode electric chillers (DMEC)

The number of DMECs of type m active at time-step t should be less than the total number of DMECs currently installed in the system.

S t,m ≤ ϕ ∑ φ=1 SD φ,m ∀m ∈ M DMEC , ∀t (4.25)
An DMEC converts electricity into either cooling power or ice. It thus has two operating modes which may be referred to as cooling mode n = COLD or ice-producing mode n = ICE. The total number of active DMECs in t is thus equal to the sum, over the two operating modes, of the number of DMECs operating in each mode.

Operation variables and constraints for the storage technologies

Let us now focus on modeling the energy storage technologies P STO at the operation level. Recall that M STO denotes the set of all types of storage devices: M STO = {(p, l) : p ∈ P STO , l = 1, . . . , L p }.

For the sake of simplicity, we will use two assumptions in what follows. First, we suppose that P STO ⊂ P CONT , i.e. that all storage technologies are of continuous nature. This means in particular that we assume that there is at most one storage device of type m = (p, 1) installed in the system and that the corresponding capacity available at the beginning of phase ϕ can be computed using the discrete design variables as Size

Step

m ∑ ϕ φ=1 SD φ,m . Second, we assume that |N in p | = |N out p | = 1, i.e.
that there is a single form of energy put into and released from the storage device. Let us denote by n in p (resp. n out p ) the form of the energy put into (resp. released from) the storage device. Note that the energy is stored within the device under a form n = n in p . Thus, for instance, for a cold storage tank (p = STO_COLD), the energy is put into and stored in the tank under the form n in p = ICE and released from the tank under the form n out p = COLD. As for a heat storage tank (p = STO_HEAT), the energy is put into, stored in and released from the tank under the form n in p = n out p = HEAT.

Variables for energy storage commodities

We first introduce two set of continuous variables: P in t,n in p ,m (resp. P out t,n out p ,m ) gives the total amount of energy of form n in p (resp. n out p ) put into (resp. released from) by the storage device of type m in t.

Moreover, STO t,m represents the amount of energy in the form n in p stored in the storage device at the beginning of time-step t.

Constraints for energy storage commodities

For any type storage device of type m, the amount of energy stored in the device at the beginning of time-step t should stay below the current installed storage capacity.

STO t,m ≤ Size

Step

m ϕ ∑ φ=1 SD φ,m ∀m ∈ M STO , ∀t (4.29)
Moreover, the amount of energy released from the storage device during t should not exceed the amount stored in the device at the beginning of t.

P out t,n out p ,m ≤ STO t,m ∀m ∈ M STO , ∀t (4.30)
The time evolution of the energy stored in the device is described through inventory balance equations. We thus have: In practice, the entering inventory of a given day in the scheduling horizon is imposed by the leaving inventory of the previous day. In our case, we do not consider each individual day of the scheduling horizon but rather a number of representative days which will not necessarily occur successively in practice. We thus impose that the entering inventory of a representative day (ϕ, d) ∈ D, STO ϕ,d,0,m , should be equal to the leaving inventory of the same selected day STO ϕ,d,H,m . This might be understood as the fact that the representative day d will be cyclically repeated w ϕ,d times in the simplified scheduling horizon used in our optimization problem for phase ϕ.

STO ϕ,

Operation constraints linking commodities and technologies

On one hand, Subsection 4.3.3 presented the operation constraints relative to the commodities using variables Inflow t,c,p or Outflow t,c,p . Each of these variables represents the aggregate amount of commodity c produced/consumed in t by all the devices belonging to technology p. Note that these Inflow or Outflow variables do not allow to follow the exact production/consumption of c by each type l = 1, . . . , L p of devices belonging to technology p.

On the other hand, Subsections 4.3.4 and 4.3.5 presented the operation constraints relative to the technologies using variables P in t,n,m and P out t,n,m . Each of these variables represents the amount of energy of form n consumed/produced in t by the devices of type m. Note that these P in and P out variables do not allow to follow the exact consumption/production by the devices of type m of the various commodities c = (n, s), s = 1, . . . , S n , corresponding to the energy form n.

In order to link the variables relative to commodities to the ones relative to technologies in our formulation, we assume that a given technology p has a single point of entry (resp. a single point of exit) for each form of energy n ∈ N in p (resp. n ∈ N out p ). The point of entry (resp. of exit) corresponding to the energy form n is connected to each commodity belonging to C in n,p (resp. C out n,p ). It is also connected to each set of devices m = (p, l), l = 1, . . . , L p . Moreover, we assume that the flows of the various commodities entering (resp. exiting) technology p at a given point are merged and shared amongst the various types of device belonging to M p .

Under these assumptions, the constraints linking commodities and technologies at the operation level are flow balance equations expressed at each entry/exit point of each technology. In the RES diagram provided in the upper part of the figure, technology p 1 consumes three commodities c 1 = (n 1 , 1), c 2 = (n 1 , 2) and c 3 = (n 2 , 1) , and produces commodity c 4 = (n 3 , 1). We thus have C in n 1 ,p 1 = {c 1 , c 2 }, C in n 2 ,p 1 = {c 3 } and C out n 3 ,p 1 = {c 4 }. The lower part of the figure displays the flows of energy, together with the corresponding operation variables, at each entry/exit point of technology p 1 . Observe how the entering flows of commodities corresponding to the same form n are merged at the corresponding entry point and how the resulting total amount of energy of form n is allocated between the various types of device. Similarly, the energy of form n produced by all types of devices is merged at the corresponding exit point before being allocated between the exiting flows of commodities of this form.

1 c3 c2 c1 c4 l 1 l 2 p 1 n 1 n 2 n 3 n 1 n 2 n 3 n 1 n 2 n 3 Outflow t,

Objective Function

Finally, while optimizing the design of an LMES, we seek to minimize the net present cost of the system which comprises the actualized design and operation costs. The objective function is thus given by: min d,h,c Sold ϕ,d,h,c (4.35)

Φ ∑ ϕ=1 α ϕ ∑ m∈M (FC m + MC ϕ,m )SD ϕ,m + ∑ c∈C R,Co SC ϕ,c Cdis c,ϕ + β ϕ ∑ d∈D ϕ w ϕ,d H-1 ∑ h=0 ∑ c∈C R EP ϕ,d,h,c Conso ϕ,d,h,c -∑ c∈C S,Sell EP ϕ,
The cost computed for each phase ϕ involves two main parts. The first part corresponds to the direct cost of the design decisions, i.e. to the cost of installing and maintaining the newly added conversion and storage devices, plus the cost of the contract negotiated with the supplier of each resource commodity in C R,Co . The second part can be seen as an estimation of the cost that will have to be paid at phase ϕ to operate the designed system so as to satisfy the energy demand at all time. This operation cost is estimated by computing the weighed sum, over all representative days in D ϕ , of the cost of the best operation schedule found for each representative day. Note that the operation cost at day (ϕ, d), is equal to the sum, over all timesteps h, of the cost of buying the resource commodities consumed during h minus the revenue obtained by selling the supply commodities during h.

Formulation of a generic technology

In Subsections 4.3.2, 4.3.4 and 4.3.5, we described how the design problem may be formulated as an MILP for a subset of energy conversion and storage technologies, focusing mostly on the ones encountered in our case studies. However, we did not cover all technologies that may be encountered in an LMES. Therefore, as a first step towards developing a tool able to deal with a larger class of LMESs, we define in this section a framework for formulating the problem for a generic technology. This framework can be seen as a template that should be customized to add new technologies within the optimization tool when necessary. Using the terminology commonly used in UML, the template technology can be considered as an abstract class and each technology mentioned above is a class which extends the abstract class.

For a given type of device m belonging to technology p, the mandatory decision variables for a general technology are the discrete design variables SD ϕ,m and the continuous variables P in t,n,m and P out t,n,m used to follow the commodity flows between the set of devices of type m ∈ M p and the other components of the LMES.

However, these variables may not be sufficient to fully describe the detailed internal functioning of a technology. We thus may have to introduce additional variables such as the number of active devices S t,m , the number of devices operating in a given mode Mode t,n,m or the amount of energy stored in the devices STO t,m . These additional auxiliary variables will be useful to model the internal behaviour of a the technology (rather than its exchanges with the other commodities). Let V t,m,a , a ∈ {1, • • • , A}, denote the set of auxiliary (integer or continuous) variables needed to model the internal functioning of a set of devices m ∈ M p at t.

The set of constraints relative to any technology can be roughly classified into three subsets.

First, the number of devices/number of discrete capacity units of a type of device m ∈ M p installed in the system should not exceed the maximum allowed number.

Φ ∑ φ=1 SD φ,m ≤ MaxNum m (4.36)
Second, in general, the maximum output power of a type of device m producing energy of type n is constrained by the number of installed devices/discrete capacity units. We cannot give a unique and explicit expression for the constraints describing how the design decisions affect the output power of a technology as different technologies may have different rules for this maximum output power: see e.g. Constraints (4.9)-(4.12) for CCHP units and Constraints (4.29)-(4.30) for storage devices. To represent these linear constraints in a generic way, we introduce a linear function O m . These constraints are formulated in a general form as follows:

O m ( Φ ∑ φ=1 SD φ,m , P out t,n,m , V t,n,m,1 , • • • , V t,n,m,A ) ≤ 0 (4.37)
Third, there is a last set of constraints to handle the relationship between the input power and output power of a type of technology. This relationship may again take many different forms. For instance, for a heat pump or a CCHP unit, the efficiency coefficients are constant and these constraints are relatively easy to formulate. For single-mode and dual-mode electric chillers, we have a set of convex non-linear performance curves for which a piece-wise linear approximation may be built. For thermal storage, the output power is directly determined by the input power and can be chosen within a range defined by the variable of the amount of stored energy. To represent these linear constraints in a generic way, we introduce a linear function E m . These constraints may also involve auxiliary variables and have the following general form.

E m (P in t,n,m , P out t,n,m , V t,n,m,1 , • • • , V t,n,m,A ) ≤ 0 (4.38)

Computational experiments

We would like to assess the computational performance of a mathematical programming solver using the MILP formulation provided in Section 4.3 at solving the LMES optimal design problem. To this aim, we carried out some numerical experiments using instances based on the three cases studies reported in Section 3.5.

Instances District cooling system in City A

The first set of instances are based on a DCS to be built in City A, located in the south of China, using as a starting point the RES diagram provided in Figure 3.3. We thus have the set of available technologies P = {SMEC, DMEC, STO_COLD} and the set of energy forms N = {ELEC, COLD, ICE}. The set of resource commodities is C R = {(ELEC, 1)}, the set of intermediate commodities C I = {(ICE, 1)} and the set of supply commodities C S = {(COLD, 1)}.

The expected lifetime of the DCS is Y = 30 years. According to the construction plan of the district, the population of the district is anticipated to increase during the first three years and stay stable afterwards. We thus introduce Φ = 3 investment phases. Phase 1 corresponds to the first year (i.e. y 1 = y 1 = 1), Phase 2 to the second year (i.e. y 2 = y 2 = 2) and Phase 3 to the last 28 years (i.e. y 3 = 3 and y 3 = 30). The annual discount rate used in the objective function to actualize the future costs is set to γ = 8%.

For each phase ϕ, we have a time series representing the predicted hourly evolution of the demand for the supply commodity (COLD, 1) throughout the year. These time series are obtained by combining historical data on the weather conditions and the typical cooling consumption in the area with forecasts on the number of buildings who will connect to the DCS. In order to highlight the upward trend of the resulting cooling demand, the total annual cooling demand and the maximum hourly demand expected for each of the three phases are provided in Table 4 As for the price of the resource commodity (ELEC, 1), it depends on the time of the day and of the commodity produced but does not vary from one day to the next. Figure 4.6 displays the evolution of the electricity price over the course of a day. It shows among others that, within a day, three types of periods may be identified. Valley periods are usually at night and early morning and correspond to a low electricity price. Peak periods happen when the demand for electricity is high, i.e. at noon and in the evening, and correspond to a high electricity price. Finally, the other periods are referred to as flat periods. It is possible to take advantage of these price variations to reduce the total energy cost of the DCS. Thus, in practice, the DCS operators often produce ice at night during valley periods, store it for a few hours and melt it to obtain chilled water at daytime during peak periods. In City A, this strategy is in fact encouraged by the grid operator which offers during valley periods an electricity price for producing ice power lower than the one for producing cooling power. In this specific case study, the price of resource commodity (ELEC, 1) thus depends on the time of the day and of the type of commodity produced.

As the price of the resource commodity does not vary with the day in the year, we select the set of representative days D ϕ using only the time series representing the cooling demand. Moreover, note that City A is located in a tropical region so that customers need cooling power during most of the year. Accordingly, the number of zero-demand days per year is relatively low: over the 3 considered phases, there is on average only 35 zero-demand days per year. The set D ϕ should thus comprise a relatively large number of days to represent as best as possible the remaining nonzero demand days. We considered the following values of |D ϕ | in our numerical experiments: {6, 14, 22, 30, 38, 58, 78}. D ϕ comprises |D ϕ | -4 typical days selected through the approach described in Subsection 4.2.2 and 4 extreme days (the day with the highest total demand, the one with the highest hourly demand, the one with the lowest non-zero total demand and the one with the lowest non-zero hourly demand). In Finally, for the resource commodity (ELEC, 1), there is a contract, to be negotiated at each phase ϕ with an outside electricity provider, fixing a limit on the maximum instantaneous power that may be bought from the grid, i.e. (ELEC, 1) ∈ C R,Co . The maximum contracted power should stay below MaxC (ELEC,1) = 30MW. We use a discrete unit of power equal to C step c = 3000kW so that MaxUnits (ELEC,1) = 10. Moreover, the price per contracted kW for is equal to 280CNY for all phase ϕ so that SC (ELEC,1),ϕ = 280 × C step c = 840kCNY. Regarding the energy conversion technologies P CONV = {SMEC, DMEC}, we have L SMEC = 3 types of single-mode electric chillers and L DMEC = 2 types of dualmode electric chillers available. For each type of device m, the installation cost FC m and the maximum production capacity P max n,m , for each n ∈ P out p , are shown in Table 4.4. The installation cost FC m of each type is assessed according to the typical cost of chillers in the market and the regional conditions for construction. Note that the chillers of type (SMEC, 1) and (SMEC,2) have the same cooling capacity P max COLD,m , but chillers of type (SMEC, 2) are less energy-efficient and less expensive than chillers of type (SMEC, 1). For each type of chiller and each commodity, the minimum output power P min n,m is equal to 0.1P For each device type m, the maximum number of devices that may be installed, MaxNum m , is equal to 10 and the maintenance cost MC ϕ,m is assumed to be equal to 0 for all phase ϕ. is assumed to be equal to 0 for all phase ϕ.

District cooling system in City B

The second set of instances are based on a DCS to be built in City B, located in the east of China, using as a starting point the RES diagram provided in The expected lifetime of the DCS is Y = 30 years and we anticipate the population of the district to increase during the first five years before getting stable. We thus consider Φ = 5 investment phases, with the first four phases lasting one year and the fifth one lasting 26 years (i.e. y 5 = 5 and y 5 = 30). The annual discount rate used in the objective function to actualize the future costs is set to γ = 8%.

Similar to City A project, we have, for each phase ϕ, a time series representing the predicted hourly evolution of the demand for the supply commodity (COLD, 1) throughout the year. In order to highlight the upward trend of the resulting cooling demand, the total annual cooling demand and the maximum hourly demand expected for each of the five phases are provided in The price of the resource commodity (ELEC, 1) features some daily periodic variations but is otherwise stationary, i.e. does not vary from one day to the next. Moreover, contrary to City A project, there is no discounted electricity price to produce ice during valley periods: the electricity price thus does not depend on the produced commodity in City B project. Figure 4.10 displays the evolution of the electricity price over the course of a day.

As the price of the resource commodity does not vary with the day in the year, we select the set of representative days D ϕ using only the time series representing the cooling demand. City B is located in a region with a temperate climate. Thus, in contrast with City A where cooling power is needed year round, the cooling demand in City B is mainly concentrated around the summer and the average number of zero-demand days per year in City B project is high: over the 5 considered phases, there is on average 204 zero-demand days per year. This is why we use a rather small number of representative days to represent in our optimization model the remaining non-zero demand days and set the value of |D ϕ | to {6, 8, 10, 12, 14}. D ϕ comprises Finally, the resource commodity (ELEC, 1) belongs to C R,Co . The maximum contracted power should stay below MaxC (ELEC,1) = 30MW. We use a discrete unit of power equal to C step c = 3000kW so that MaxUnits (ELEC,1) = 10. Moreover, the price per contracted kW for is equal to 440CNY for all phase ϕ so that SC (ELEC,1),ϕ = 440 × C step c = 1.32MCNY. Regarding the energy conversion technologies P CONV = {SMEC, DMEC}, there are L STDC = 2 types of single-mode electric chillers and L ICEC = 1 type of dual-mode electric chiller that may be installed. These types correspond to the large-capacity chillers (SMEC, 1), (SMEC, 2) and (DMEC, 1) used for City A project: the corresponding performance curves at an ambient temperature of 30 For each device type m, the maximum number of devices that may be installed, MaxNum m , is equal to 10 and the maintenance cost MC ϕ,m is assumed to be equal to 0 for all phase ϕ. Finally, we have a single energy storage technology and this one is of continuous nature: P STO = {STO_COLD} ⊂ P CONT . The capacity of the thermal storage tank m = (STO_COLD, 1) to be installed may namely take any value between 0 MWh and MaxSize m = 180MWh. We use a discrete capacity unit Size step m equal to 9kWh so that MaxNum m = 20 . The cost for installing a kWh of ice storage capacity is 130CNY so that FC m = 130 × Size step m = 1.17MCNY . The maintenance cost MC ϕ,m is assumed to be equal to 0 for all phase ϕ.

Trigeneration system in City C

The last set of instances are related to the trigeneration system to be built in City C, using as a starting point the RES diagram provided in Figure 3.4. We thus have the set of available conversion technologies P CONV = {CCHP, BOILER, SMEC, DMEC, ASHP}, the set of available storage technologies P STO = {STO_COLD, STO_HEAT} and the set of energy forms N = {ELEC, GAS, HEAT, COLD, ICE}. The set of resource commodities is C R = {(ELEC, 1), (GAS, 1)}, the set of intermediate commodities C I = {(ELEC, 2), (HEAT, 1), (ICE, 1)} and the set of supply commodities C S = {(HEAT, 2), (COLD, 1)}.

The expected lifetime of the trigeneration system is Y = 30 years and the population of the district is expected to increase in the first four years before getting stable. Therefore, we consider Φ = 4 construction phases, with the first three phases lasting one year and the fourth one lasting 27 years (i.e. y 4 = 4 and ȳ4 = 30). The annual discount rate of this project is set to γ = 8%.

For each phase ϕ and each supply commodity, we have a time series describing the predicted hourly evolution over the year. The total yearly demand and the maximum hourly demand are given in Table 4.7 for (HEAT, 2) and (COLD, 1): in both cases, we observe a clear upward trend over the years. The heat dissipation is integrated into the numerical values of energy demand for (HEAT, 2). Moreover, City C is located in central China where the climate is temperate. The cooling demand and the heating demand of the neighboring residents thus display a seasonal character and occur in two non-overlapping periods during a year. More precisely, the cooling supply period covers six months, from the 1st May to the 31st October and the heating supply period covers four months, from the 15th November to the 15th March of the next year.

The price of the resource commodity (ELEC, 1) bought from the national grid varies within the day. Four periods may be distinguished: valley, flat, high and peak periods. The valley period is from 23:00 to 7:00 the next day and corresponds to a unit price of 0.27 CNY/kWh. The high period is from 9:00 to 16:00 with a price 0.85 CNY/kWh. The peak period lasts for two hours from 20:00 to 22:00 with a price 1.03 CNY/kWh. The rest of the day is termed flat period and corresponds to a unit price of 0.57 CNY/kWh. The variation of the electricity price is shown in Figure 4.12. The buying price of the resource commodity (GAS, 1) natural gas is assumed to stay stable during the lifetime of the system unit and is set to is 0.27 CNY/kWh. As the price of the resource commodities does not vary with the day in the year, we select, for each phase ϕ, a set of representative days using only the time series corresponding to the predicted hourly evolution of the cooling demand and of the heating demand throughout the year. We thus select a set of |D ϕ | ∈ {12, 16, 20, 24, 28, 32, 36} representative days: |D ϕ | -8 typical days selected using the procedure described in Subsection 4.2.2, 4 extreme days for the cooling demand and 4 extreme days for the heating demand. Each day is divided into H = 24 time-steps and each time-step is of length ∆t = 1 hour. Figure 4.13 displays the 4 extreme days selected for each type of demand at phase 4.

Finally, the maximum instantaneous electric power that may be bought from the electricity power is limited through a contract which may be renegotiated at the beginning of each phase. We thus have C R,Co = {(ELEC, 1)}. The maximum contracted power should stay below MaxC (ELEC,1) = 30MW. We use a discrete unit of power equal to C step c = 0.5MW so that MaxUnits (ELEC,1) = 60. Moreover, the price both modes is 30% so that P min COLD,(DMEC,1) = 1.1MW and P min ICE,(DMEC,1) = 0.72MW. Its conversion efficiency is load-and temperature dependant. It is depicted by two sets of convex performance curves providing the consumed electric power as a function of the produced cooling power (or produced ice) under different ambient temperatures. The performance curves of both modes are shown in Figure 4.14. The purchase price for one unit of (DMEC, 1) is 1.8 MCNY and the total installation cost is FC (DMEC,1) = 4.4MCNY. There is one type of air-source heat pump available to in the catalog: L ASHP = 1. Its heating and cooling capacities are equal to P max COLD,(ASHP,1) = P max HEAT,(ASHP,1) = 1.6MW. There is no mininum loading rate so that P min COLD,(ASHP,1) = P min HEAT,(ASHP,1) = 0MW. The conversion efficiency of the heating mode depends on the ambient temperature: we have η t,HEAT,(ASHP,1) = 0.064T outdoor,t + 1.75, in which T outdoor,t is the predicted ambient temperature at time-step t. The cooling mode conversion efficiency is constant: η t,COLD,(ASHP,1) = 2 for all t. The purchase price for a unit of (ASHP, 1) is 1.8MCNY and the total installation cost FC (ASHP,1) = 4.4 MCNY.

The description of the available energy conversion devices is summarized in Table 4.8. The installation cost FC m comprises both the purchase cost PC m and the construction cost CC m of a device. The annual maintenance cost of any installed device is set to 2.5% of its purchase cost for each year of the lifetime of the system. This gives MC m,ϕ = 0.025

Results

In this subsection, we report the results of the computational experiments carried out on the instances based on the three case studies described above.

All tests were carried out on a machine with an Intel Xeon 2.90GHz processor with 8 cores and 16GB RAM, running under Windows 7. All the algorithms were implemented in Python. The direct resolution of the LMES design problem was done by the mathematical programming solver CPLEX 12.8 using the MILP formulation introduced in Section 4.3. We set the solver to use up to 8 threads and a maximum RAM of 11.2GB (70% of the machine's RAM) during the resolution, impose a time limit of 7200s and otherwise use the solver default settings.

The numerical results are provided in Table 4.9, Table 4.10 and Table 4.11 for the instances corresponding to City A DCS project, City B DCS project, and City C trigeneration project.

For each instance, we display the following information:

• #Var: the total number of continuous or integer decision variables involved in the MILP formulation.

• #IntVar: the total number of integer decision variables.

• #Cons: the total number of constraints.

• Gap LP : the integrality gap, i.e. the percentage relative difference between the lower bound provided by the linear relaxation of the problem and the optimal solution of the problem. Note that the exact optimal objective value is found using the hierarchical decomposition algorithm and/or the generalized benders decomposition algorithm, which will be introduced in Chapters 5 and 6.

• Time: the total computation time of the algorithm before a guaranteed optimal solution of the design problem is found or the time limit is reached.

• Gap MIP : the percentage residual gap between the best lower bound and the best integer feasible solution of the design problem found by the algorithm within the time limit. Note that Gap = 0.00% in case a guaranteed optimal solution of the problem was found within the time limit.

Results from Tables 4.9-4.11 first show that the direct resolution of design problems by a mathematical solver never provides a guaranteed optimal solution within the computation time limit but that near-optimal solutions (i.e. solutions displaying a residual gap below 1.0%) can be found for most instances. However, numerical difficulties arise when the number of representative days |D ϕ | and consequently the MILP size increase. This can be seen for the City A project instance corresponding to |D ϕ | = 78 for which the residual gap is large (i.e. equal to 2.83%). Moreover, out-of-memory issues were encountered during the resolution for the City B project instances corresponding to |D ϕ | ∈ {12, 14} and for the City C project instances corresponding to |D ϕ | ∈ {12, 16, 24, 28, 32, 36}.

Furthermore, we found that the system deployment plan given by the best found feasible solution may significantly change with the number of representative days considered for each phase in the instance. This may be due to the fact that the number of representative days is not large enough to correctly assess the operation cost (so that the design decisions depend heavily on the subset of selected representative days) or to the fact that the best found feasible solution is not the optimal one due to the existence of non-zero residual gap Gap MIP . In any case, we have doubts about the quality of the system deployment plans obtained through this direct resolution of the MILP formulation by a solver.

|D

Conclusion

This section investigated how the problem of optimally designing an LMES over a multi-phase investment horizon may be modeled and formulated as an MILP.

We first presented the modeling assumptions relative to the definition of the investment phases and the use of representative days to estimate the operation cost. We also explain how the convexity of the performance curves of the technologies displaying a part-load efficiency may be exploited to build operation schedules at an aggregate level (instead of operation schedules for each individual device). The main objective of this first section was to drastically reduce the number of variables and constraints needed to formulate the problem as a mathematical program.

Thanks to this, we were able to propose an MILP formulation of computationally tractable size for the problem. Note that the proposed MILP formulation relies on some additional assumptions such as: all conversion (resp. storage) technologies are of discrete (resp. continuous) nature, each technology consumes energy under a single form and all conversion technologies have linear or convex energy conversion performance curves. The proposed MILP formulation enabled us to model the three case studies reported in Chapter 3 and could be applied to a large number of LMES design problems. However, further research might be needed in order to extend this MILP formulation so as to be able to design a broader class of LMESs. In particular, it would be interesting to study LMESs involving renewable energy sources, conversion technologies consuming energy under two different forms (such as absorption chillers), conversion technologies with non-convex performance curves or inter-seasonal thermal storage in which energy is stored for long periods spanning several months. This is however left for future work.

Finally, we assessed the numerical efficiency of a mathematical solver at providing optimal solutions by directly solving the problem using the proposed MILP formulation. Our computational experiments carried out on a set of 19 large-size instances coming from the three reported real-life case studies indicated many numerical difficulties, in particular out-of-memory issues for the largest instances. This motivated us to study the development of more advanced solution approaches using decomposition-based algorithms which will be described in the next two chapters.

Chapter 5

Hierarchical decomposition algorithm 5.1 Introduction

The problem modeling proposed in Chapter 4 leads to the formulation of a large size MILP which may be solved directly by a commercial MILP solver. The computational results obtained with this direct resolution approach for our three case studies are displayed in Tables 4.11-4.11. They show that, even with the smallest number of representative days, the resolution process cannot converge before the time limit, and that, for some instances, the resolution is terminated due to out-ofmemory issues. These numerical difficulties can mainly be attributed to the large number of integer variables and constraints introduced in the formulation to build the operation schedule for each representative day.

In the present chapter, we discuss a decomposition algorithm able to solve large instances of the LMES optimal design problem. This algorithm exploits the fact that the MILP problem introduced in Chapter 4 shows a special structure. Namely, the variables of the MILP model defined in Chapter 4 can be classified into two classes. The design variables refer to the decisions related to the sizing of technologies and to the contracts determining the maximum allowed power consumption of resource commodities, which are made at the beginning of each phase. The operation variables refer to decisions on the input-output of devices, inflow-outflow of commodities, the operating status of devices in each time-step. The decomposition algorithm investigated here is based on the key observation that, once the values of the design variables are determined for all construction phases, the original problem can be decomposed into a series of independent sub-problems. Each of these sub-problems takes the current value of the design variables as input and seeks to minimize, for the corresponding system layout, the operation cost for one representative or extreme day. The solution approach investigated in this chapter exploits this hierarchical structure and is based on a bi-level decomposition of the problem. The original problem is thus decomposed into an upper-level design problem and a set of lowerlevel operation scheduling sub-problems of smaller size. The algorithm proposed in this section can be seen as an extension of the hierarchical decomposition algorithm introduced by [START_REF] Yokoyama | Optimization of energy supply systems by MILP branch and bound method in consideration of hierarchical relationship between design and operation[END_REF]. The proposed extended algorithm mainly aims to be more computationally efficient at solving problems involving multiple deployment phases of the LMES and short-term thermal storage capacity.

This chapter is organized as follows. In Section 5.2, we formulate the original problem in a concise and compact way to illustrate the bi-level hierarchical structure of the problem. In Section 5.3, the original hierarchical decomposition algorithm is presented. In Section 5.4, we investigate the proposed extension of the original algorithm. The computational results obtained for our three case studies are shown in Section 5.5.

Compact reformulation and hierarchical structure of the problem

In order to explain the hierarchical decomposition algorithm in an easier and more concise way, we introduce in the present section a compact formulation of our problem.

In Chapter 4, we introduced two sets of design variables: integer variables SD φ,m to represent the number of devices or the number of discrete capacity units of type m ∈ M added in the system at the beginning of phase ϕ, and integer variables Cdis ϕ,c to represent the number of discrete units of power contracted for the resource commodity c ∈ C R,Co in phase ϕ. In order to obtain a more compact MILP formulation, we form a single integer variable vector by juxtaposing these design variables in a certain order. Therefore, we order the elements in the set M from m 1 to m |M| and the element in the set C R,Cp from c

1 to c |C R,Co | . Let ν ϕ = |M| + |C R,Co |.
The system layout at phase ϕ is then represented in the compact formulation by the following integer vector:

x ϕ = ( ϕ ∑ φ=1 SD φ,m 1 , • • • , ϕ ∑ φ=1 SD φ,m |M| , Cdis ϕ,c 1 , • • • , Cdis ϕ,c |C cont | ) ∈ Z ν ϕ (5.1)
x ϕ gives the total number of devices / discrete capacity units of each type of device in M present in the system at the beginning of phase ϕ and the number of basic instantaneous power units contracted for each resource commodity in C R,Co at phase ϕ. The integer variable vector x = (x 1 , ..., x Φ ) ∈ Z ν , with ν = ∑ Φ φ ν φ , describes the deployment plan over the whole investment horizon.

We denote by D = {(ϕ, d)|ϕ ∈ {1, ..., Φ}, d ∈ D ϕ } the set of all representative days. For each day (ϕ, d) ∈ D, we introduce the vector y ϕ,d to represent in a compact way the continuous operation variables needed to describe the operation schedule during this day. y ϕ,d thus stands for all the Inflow, Outflow, Conso, Sold, P in , P out and Sto variables relative to day (ϕ, d). Similarly, we introduce the vector z ϕ,d to represent in a compact way the integer operation variables needed to describe the operation schedule during this day. z ϕ,d thus stands for all the S and Mode variables relative to day (ϕ, d). Let µ (resp. λ) denote the dimension of vector y ϕ,d (resp. z ϕ,d ).

Note that, in Chapter 4, we used bold letters for the decision variables and italic letters for the input parameters of the MILP formulation. Here, for the sake of readability, we change our convention. We thus use a bold letter to denote a vector of variables and an italic letter to denote a single variable. Input parameters do not appear frequently in the compact formulation, but when they appear, we will give detailed explanation to avoid ambiguity.

Let us now have a closer look at the constraints introduced in the MILP formulation described in Chapter 4. We first have a set of design constraints, i.e. constraints involving only design decisions: see Constraints (4.1)-(4.2). We then have a set of constraints coupling together the design decisions and the operation decisions relative to a given representative day (ϕ, d): see e.g. Constraints (4.9), (4.15), (4.19), (4.22), (4.25) for the energy conversion technologies, Constraints (4.29) for the energy storage technologies and Constraints (4.5) for the resource commodities in C R,Co . Finally, all the remaining constraints may be seen as operation constraints involving only operation decisions relative to a given representative day (ϕ, d).

As a consequence, the LMES design problem can be reformulated compactly as follows:

min Z CP = f 0 (x) + ∑ (ϕ,d)∈D f ϕ,d (y ϕ,d , z ϕ,d ) (5.2) h(x) ≤ 0 (5.3) g ϕ,d (x ϕ , y ϕ,d , z ϕ,d ) ≤ 0 ∀(ϕ, d) ∈ D (5.4) l ϕ,d (y ϕ,d , z ϕ,d ) ≤ 0 ∀(ϕ, d) ∈ D (5.5)
x ∈ Z ν (5.6)

y ϕ,d ∈ R µ ∀(ϕ, d) ∈ D (5.7) z ϕ,d ∈ Z λ ∀(ϕ, d) ∈ D (5.8)
The objective function (5.2) involves two terms.The first term f 0 (x) corresponds to the total net present design cost of the deployment plan x. The second term computes the sum, over all selected days, of the net present operation cost f ϕ, d (y ϕ,d ,z ϕ,d ) of the schedule computed for day (ϕ, d) and described by (y ϕ,d , z ϕ,d ). Constraints (5.3) represent the design constraints. Constraints (5.4) represent the coupling constraints which link together design and operation variables. Finally, Constraints (5.5) stand for the operation constraints which are expressed using only operation variables. Constraints (5.6)-(5.8) give the domain definition of the variables. Problem (5.2)-(5.8) will be referred to as the Complete Problem (CP) in what follows.

The hierarchical decomposition algorithm uses as the upper level problem a relaxation of CP in which the design variables are kept integer but the discrete operation variables z ϕ,d are relaxed to continuous variables. The resulting problem is referred to as the semi-relaxed problem (SRP). SRP has the same number of variables and constraints as CP but a significantly smaller number of integer variables. Its resolution by an MILP solver should thus be easier than the one of CP.

min Z SRP = f 0 (x) + ∑ (ϕ,d)∈D f ϕ,d ( ỹϕ,d , zϕ,d ) (5.9)
h(x) ≤ 0 (5.10)

g ϕ,d (x ϕ , ỹϕ,d , zϕ,d ) ≤ 0 ∀(ϕ, d) ∈ D (5.11) l ϕ,d ( ỹϕ,d , zϕ,d ) ≤ 0 ∀(ϕ, d) ∈ D (5.12) x ∈ Z ν (5.13) ỹϕ,d ∈ R µ ∀(ϕ, d) ∈ D (5.14) zϕ,d ∈ R λ ∀(ϕ, d) ∈ D (5.15)
The hierarchical decomposition algorithm is based on the key observation we mentioned above: when the system deployment plan is determined, i.e., when the value of vector x is fixed, the problem can be decomposed into a series of independent lower-level scheduling sub-problems. Let x ♯ denote a fixed system deployment plan. The operation sub-problem relative to representative day (ϕ, d) ∈ D, denoted by OP ϕ,d (x ♯ ϕ ), is defined as follows:

min Z ϕ,d (x ♯ ϕ ) = f ϕ,d (y ϕ,d , z ϕ,d ) (5.16) g ϕ,d (x ♯ ϕ , y ϕ,d , z ϕ,d ) ≤ 0 (5.17) l ϕ,d (y ϕ,d , z ϕ,d ) ≤ 0 (5.18) y ϕ,d ∈ R µ (5.19) z ϕ,d ∈ Z λ (5.20) OP ϕ,d (x ♯ ϕ )
is an MILP problem involving a number of operation variables and constraints much smaller that the one involved in CP. It is thus much easier to solve than the original complete problem CP.

Decomposition algorithm

Outline of the algorithm

The hierarchical decomposition algorithm proposed by [START_REF] Yokoyama | Optimization of energy supply systems by MILP branch and bound method in consideration of hierarchical relationship between design and operation[END_REF] exploits the hierarchical structure of the problem discussed in Section 5.2. We note that a prerequisite to apply this algorithm is that all design variables should be discrete. This restriction comes from the customized branch-and-bound algorithm used to solve the upperlevel design problem, which will be explained in detail in what follows.

At the upper level of the algorithm, SRP is solved by a branch-and-cut (B&C) algorithm. Since all the design variables remain discrete in SRP, each integer feasible solution (x ♯ , ỹ♯ , z♯ ) of SRP found during the B&C tree search corresponds to a potential deployment plan described by the integer vector x ♯ and complying with the design constraints (5.3) and (5.6). However, as the integer scheduling variables z ϕ,d have been relaxed to continuous variables zϕ,d in SRP, an integer feasible solution (x ♯ , ỹ♯ , z♯ ) of SRP is not necessarily an integer feasible solution of CP. Moreover, the feasible region of CP is a subset of the feasible region of SRP. Therefore, the objective value of SRP,

Z ♯ SRP = f 0 (x ♯ ) + ∑ (ϕ,d)∈D f ϕ,d ( ỹ♯ ϕ,d , z♯ ϕ,d
), provides a lower bound to the actual net present cost Z ♯ CP of the deployment plan x ♯ . Hence, during the B&C tree search, each time an integer feasible solution (x ♯ , ỹ♯ , z♯ ) of SRP is found, we should check its feasibility for CP with respect to the operation constraints (5.4)-(5.5), (5.7)-(5.8) and compute its accurate net present cost Z ♯ CP . This can be done by solving a sequence of independent lower level operation sub-problems. To be more precise, we solve, for each day (ϕ, d) ∈ D, the problem OP ϕ,d (x ♯ ϕ ), in which the operation decision variables z ϕ,d ∈ Z λ are integer. This allows us to first check the feasibility of the system layout x ♯ ϕ for each selected day. Moreover, if OP ϕ,d (x ♯ ϕ ) is a feasible problem, we solve it to optimality to obtain an optimal operational schedule (y

♯ ϕ,d , z ♯ ϕ,d ), together with its corresponding exact operation cost f ϕ,d (y ♯ ϕ,d , z ♯ ϕ,d ). If all sub-problems OP ϕ,d (x ♯ ϕ ), (ϕ, d) ∈ D are feasible, (x ♯ , y ♯ , z ♯ ) is a feasible solu- tion of the original CP with a net present cost Z ♯ CP = f 0 (x ♯ ) + ∑ (ϕ,d)∈D f ϕ,d (y ♯ ϕ,d , z ♯ ϕ,d
(x) + ∑ (ϕ,d)∈D f ϕ,d ( ỹϕ,d , zϕ,d ) ≤ Z best
CP to the formulation of SRP.

Note that, whether the feasible solution (x ♯ , y ♯ , z ♯ ) is accepted as a new incumbent of CP or not, the potential deployment plan (x ♯ , ỹ♯ , z♯ ) is rejected in the B&C search tree used to solve SRP at the upper level. Namely, the solver should not take the integer feasible solution (x ♯ , ỹ♯ , z♯ ) as the new incumbent of SRP but rather should continue branching from this node.

This rejection is necessary for two reasons. First, if the current integer feasible solution (x ♯ , ỹ♯ , z♯ ) is not rejected as an incumbent of SRP, the upper-level B&C algorithm may use the corresponding objective value Z ♯ SRP as the current best known upper bound and decide to prune nodes in the subsequent iterations based on this value. However, Z ♯ SRP is a lower bound of the actual net present cost and corresponds only to an approximate estimation of the cost performance of the potential deployment plan x ♯ . It therefore should not be used to compare x ♯ with other potential deployment plans and to decide which nodes to prune in the B&C tree. Second, let n ♯ denote the node of the B&C tree on which the potential deployment plan (x ♯ , ỹ♯ , z♯ ) is found and let F ♯ be the subset of the feasible domain of SRP defined by the branching decisions made up to node n ♯ . For any solution (x, ỹ, z) in F ♯ , we have Z SRP ≥ Z ♯ SRP . It is therefore not possible to find in F ♯ an integer feasible solution of SRP with a value Z SRP lower than Z ♯ SRP . However, as Z SRP is only an approximation of the actual net present value of a deployment plan, we may find in

F ♯ an integer feasible solution (x • , ỹ• , z• ) of SRP corresponding to an integer feasi- ble solution (x • , y • , z • ) of CP such that Z • CP < Z ♯ CP .
In other words, we should not decide to stop the exploration in F ♯ by considering only the fact that (x ♯ , y ♯ , z ♯ ) is integer feasible for CP as better integer feasible solutions of CP may be found in F ♯ . Rejecting (x ♯ , ỹ♯ , z♯ ) as the new incumbent solution for SRP thus prevents the solver from using the value Z ♯ SRP to wrongly prune nodes in the search tree, from prematurely stopping the exploration in F ♯ and thus from neglecting other integer feasible solutions of CP contained in F ♯ . This rejection is realized by a built-in function of the mathematical solver

In fact, the upper-level B&C algorithm should never accept any incumbent of SRP in the B&C search tree and the decision to prune nodes in this tree should only be made by comparing the corresponding lower bound with the value Z best CP of the incumbent of CP. More precisely, a node in the B&C search tree of the upper-level problem should be pruned only if its corresponding deployment plan (x, ỹ, z) is such that Z SRP > Z best CP or if it is infeasible with respect to SRP. Therefore, a node n with objective value Z SRP < Z best remains active and may be used for branching even if (x, ỹ, z) is an integer feasible solution of SRP. In this case, the branching will be carried out by a built-in function of the mathematical solver on an design decision variable x ϕ,i which already has an integer value at node n. This involves that an integer potential deployment plan similar to the one found at node n may be found while exploring the sub-tree of the B&C tree rooted in n. In order to avoid evaluating multiple times the same deployment plan, we add to the formulation of SRP a nogood cut ensuring that the integer potential deployment plan found at a node of the B&C search tree is excluded from the feasible domain of SRP. This will be explained in more detail in the next two sections.

When the B&B algorithm in the upper level problem terminates, i.e. when all nodes have been closed, the incumbent of CP (x best , y best , z best ) gives an optimal solution of CP.

0 ( 3 2 , 0) 1 (1, 2 3 ) 3 infeasible y ≤ 0 4 (1, 1) 5 (0, 2) 11 infeasible y ≤ 1 12 (0, 2) 13 (0, 2) y ≤ 2 14 infeasible y ≥ 3 y ≥ 2 x ≤ 0 6 (1, 1) 7 (1, 1) y ≤ 1 8 (1, 2) 9 (1, 2) y ≤ 2 10 infeasible y ≥ 3 y ≥ 2 x ≥ 1 y ≥ 1 x ≤ 1 2 infeasible x ≥ 2

Incumbent rejection

In this subsection, we use a small example to clarify what happens when the solver rejects the current integer feasible solution as incumbent and continues branching from a node where an integer feasible solution has been found. We consider the following small integer program. min y (5.21)

s.t. 4 3 x + y ≥ 2 (5.22) - 1 3 x + y ≤ 2 (5.23) x ≤ 3 2 (5.24) (x, y) ∈ Z 2 (5.25)
The feasible domain of the linear relaxation of (5.21)-(5.25) is a triangle defined by the vertices (0, 2), (1.5, 0) and (1.5, 2.5). It comprises three integer feasible solutions: (1, 1), (0, 2) and (1, 2). As the objective function consists in minimizing the value of the integer variable y, the optimal integer solution is (1, 1).

Let us now solve Problem (5.21)-(5.25) with a B&C algorithm in which no integer feasible solution is accepted as incumbent during the solving process. The corresponding search tree is shown in Figure 5.1.

We first note that, if we accept the current integer feasible solution as the incumbent when suitable, the B&C algorithm terminates at Node 4 with the optimal solution (1, 1) and optimal value 1. However, if we reject the integer feasible solution (1, 1) as incumbent, the algorithm will continue and will branch from Node 4 by adding bounds on the variable x, which already has an integer value at Node 4. Contrary to the traditional branching operation applied on fractional-valued variables, this branching will not exclude the current integer feasible solution. Therefore, we can observe that Node 6, a descendent node of Node 4, also provides the integer feasible solution (1, 1). The branching actions will stop either when the LP relaxation is infeasible (Nodes 10, 11 and 14) or when the value of each variable is restrained to take a single value (Nodes 7, 9 and 13), i.e. when the feasible domain of the problem has been partitioned into subsets involving at most one integer feasible solution. As a result, the modified B&C tree explores more nodes than a standard B&C tree and returns no feasible solution at the end of exploration. Moreover, all the integer feasible solutions contained in the feasible space are considered as integer feasible solution on at least one node of the search tree. Thanks to this small example, we see how rejecting the incumbent still guarantees the finite and optimal convergence of the B&C algorithm. First, even if we reject all found integer feasible solutions as incumbent, the B&C search tree terminates in a finite number of steps, when enough branching decisions have been made to ensure that the feasible domain has been partitioned into subsets involving at most one integer feasible solution, so that the termination of the algorithm is guaranteed theoretically. Second, all the integer feasible solutions are explored at least once during the B&C search tree. As we can see in the example, even though an integer feasible solution is found in Node 4, all the other integer feasible solutions contained in the feasible region defined by the branching decisions made up to Node 4 appear on at least one node of the sub-tree rooted at Node 4. This ensures that the optimal solution is found and evaluated at least once during the B&C search.

Let us now come back to the hierarchical decomposition algorithm mentioned above.

The simple example provided above first shows that even if we do not cut off open nodes based on a comparison between their lower bound and the current incumbent, the B&C algorithm with a systematic incumbent rejection still terminates in a finite number of steps, basically after having carried out a trivial enumeration of all integer feasible solutions. The hierarchical decomposition algorithm can be seen as an extension of this trivial enumerating algorithm in which we close nodes by comparing their lower bound with the current best known integer value for CP. The hierarchical decomposition algorithm will thus terminate in a finite number of steps but is likely to be more computationally efficient that the basic algorithm based on a full enumeration discussed above.

Moreover, thanks to the systematic incumbent rejection, when an integer feasible solution of SRP has been found at a node, we will continue branching from that node. This ensures that all the integer feasible solutions of SRP, and consequently all integer feasible solutions of CP, compatible with the branching decisions made up to that point will be explicitly considered and evaluated through the resolution of a series of operation sub-problems. The only exception corresponds to the case where the lower bound computed at a node is larger than the current best known integer value for CP. In this case, we will close the node and by doing so, we will eliminate some integer feasible solutions of CP. But as these solutions cannot be better than the current best known solution of CP, this does not affect the optimal convergence of the algorithm.

No-good cuts

However, the systematic rejection of the incumbent leads to some loss of computational efficiency as the same integer feasible solution of SRP (and the same integer feasible solution of CP) may be considered and evaluated multiple times during the resolution process.

To prevent this from happening, 'no-good cuts' may be added to the formulation of SRP [START_REF] Wolsey | Integer programming[END_REF]. No-good cuts are applied on a set of binary variables, so we need to carry out a binary expansion of the integer design variables. More precisely, for each component x ϕ,i of the vector x ϕ , we introduce a set of additional binary variables, xbin ϕ,i,θ , with θ ∈ {0, ..., Θ ϕ,i }, and define the binary expansion of x ϕ,i as:

x ϕ,i = Θ ϕ,i ∑ θ=0 2 θ xbin ϕ,i,θ (5.26)
Here, xbin ϕ,i,θ = 1 if coefficient 2 θ is used to compute the value of x ϕ,i and 0 otherwise. Note that Θ ϕ,i is set to Θ ϕ,i = ⌈log 2 x max ϕ,i ⌉, where x max ϕ,i is an upper bound of the variable x ϕ,i . This binary expansion ensures that all possible integer values of x ϕ,i within its range has a one-to-one correspondence to the binary representation

(xbin ϕ,i,0 , • • • , xbin ϕ,i,Θ ϕ,i ).
Now, if a given deployment plan x 1 = (x 1 1 , ..., x 1 Φ ) is rejected, we add a 'no-good cut' using the binary expansion variables to exclude this potential deployment plan from the feasible space of SRP. This cut is obtained by partitioning the set of all indices used in the binary expansion of each vector x ϕ , i.e., the set

I ϕ = {(i, θ), i ∈ {1, • • • , ν ϕ }, θ ∈ {0, • • • , Θ ϕ,i
}}, into two sub-sets according to the value of the binary expansion of the given deployment plan xbin 1 :

I 0 ϕ = {(i, θ) ∈ I|xbin 1 ϕ,i,θ = 0} and I 1 ϕ = {(i, θ) ∈ I|xbin 1 ϕ,i,θ = 1}.
Using this notation, we formulate the no-good cut as follows:

Φ ∑ φ=1   ∑ (i,θ)∈I 0 φ xbin φ,i,θ + ∑ (i,θ)∈I 1 φ (1 -xbin φ,i,θ )   ≥ 1 (5.27)
Constraint (5.27) enforces that any potential deployment plan considered differs from x 1 with respect to at least one of the binary variables involved in the binary expansion of x. Adding Constraint (5.27) to the formulation of SRP thus ensures that any potential deployment plan encountered during the upper level branch-and-cut search will appear only once, which allows us to avoid some repetitive computations.

Algorithmic improvement

The hierarchical decomposition algorithm discussed in Section 5.3 does not always outperform the generic B&C algorithm embedded in CPLEX 12.8 at solving the design problem, as will be shown by the numerical results reported in Section 5.5. There are two reasons for this inefficiency. Firstly, the presence of multiple construction phases leads to an increase in the number of potential deployment plans to be explored and consequently to an increase in the computational burden. Second, since all the operation discrete variables are relaxed to continuous variables in SRP, an integer feasible solution of SRP is not necessarily feasible with respect to the original problem CP. This means that the upper level algorithm explores a significant number of potential deployment plan which are feasible for SRP but infeasible for CP. In this section, we will investigate two ways to handle these two factors which slow down the resolution approach.

Leveraging previously computed results for the lower level subproblems

Part of the loss of efficiency comes from the fact that, over the course of the hierarchical decomposition algorithm, the solver often solves the same lower level subproblems repetitively. This repetition is due to the multi-phase nature of CP. Namely, even if two potential deployment plans found in the upper level B&B search tree differ from each other, they may contain the same system layout for a subset of construction phases. More precisely, two potential deployment plans x 1 = (x 1 1 , ..., x 1 Φ ) and x 2 = (x 2 1 , ..., x 2 Φ ) may be such that x 1 ̸ = x 2 and x 1 ϕ = x 2 ϕ for some but not all phases ϕ ∈ {1, ..., Φ}. Therefore, for all ϕ such that x 1 ϕ = x 2 ϕ and all days d ∈ D ϕ , we will have OP ϕ,d (x 1 ϕ ) = OP ϕ,d (x 2 ϕ ) and the computation of the exact net present cost of x 1 and x 2 , Z 1 CP and Z 2 CP , will involve solving a set of identical operation subproblems. The corresponding loss of computational efficiency is all the most important in our case that these sub-problems are formulated with a scheduling horizon comprising H > 1 time periods and are thus MILP problems of a significantly larger size than the ones used in [START_REF] Yokoyama | Optimization of energy supply systems by MILP branch and bound method in consideration of hierarchical relationship between design and operation[END_REF] to formulate the operation lower level subproblems. This comes from the fact that [START_REF] Yokoyama | Optimization of energy supply systems by MILP branch and bound method in consideration of hierarchical relationship between design and operation[END_REF] do not consider short-term (intraday) energy storage and thus formulate operation sub-problems involving only a single scheduling period (i.e. a single hour). In contrast, in our case, the existence of energy storage leads to the formulation of inventory balance equations, which couple together the H time steps of a selected day and force us to consider multi-period scheduling sub-problems.

In order to avoid this waste of computation capacity and thus reduce the computational burden of the algorithm, we propose to improve the implementation of the hierarchical decomposition algorithm.

Excluding infeasible system layouts

A first improvement deals with the case in which the system layout proposed by a given potential deployment plan x 1 = (x 1 1 , ..., x 1 Φ ) is infeasible at some phase ϕ, i.e. is such that at least one operation sub-problem OP ϕ,d (x 1 ϕ ) at phase ϕ is infeasible. In this case, we add a no-good cut to exclude all potential deployment plans x such that x ϕ = x 1 ϕ . This cut is obtained by using the same partition of the set of all indices as the one used for the no-good cuts (5.27), i.e., the sets

I 0 ϕ = {(i, θ) ∈ I|xbin 1 ϕ,i,θ = 0} and I 1 ϕ = {(i, θ) ∈ I|xbin 1 ϕ,i,θ = 1}.
The no-good cut added to the formulation of SRP to exclude the infeasible system layout x 1 ϕ is formulated as follows:

∑ (i,θ)∈I 0 ϕ xbin ϕ,i,θ + ∑ (i,θ)∈I 1 ϕ (1 -xbin ϕ,i,θ ) ≥ 1 (5.28)
Constraint (5.28) enforces that any potential system layout considered for phase ϕ differs from x 1 ϕ with respect to at least one of the binary variables involved in the binary expansion of x ϕ . Adding Constraint (5.28) to the formulation of SRP thus ensures that any potential deployment plan encountered during the upper level branch-and-cut search will be considered as a feasible solution of SRP only if at least one of the components of x ϕ has a value different to the one it has in x 1 ϕ .

Storing the optimal value of previously solved operation sub-problems

In order to avoid repetitively solving the same operation sub-problems multiple times over the course of the algorithm, we also propose to record in a hash table the optimal value of any feasible operation sub-problem solved during the algorithm. Thus, let x 1 be a potential deployment plan found during the upper level branchand-cut search. For each phase ϕ, if the sub-problem OP ϕ,d (x 1 ϕ ) relative to day (ϕ, d) is feasible, its objective value Z ϕ,d (x 1 ϕ ) is stored in the hash table with a key (x 1 ϕ , ϕ, d). This information can be latter used in the algorithm to save some computational effort. Namely, suppose that x 2 = (x 2 1 , . ) as an MILP and record its optimal value in the hash table if it is feasible. The searching time in a hash table is short and independent of the table size, it is thus negligible as compared to the time needed to solve an operation sub-problem as an MILP.

Valid inequalities

As explained in Section 5.3, an integer feasible solution of SRP (x ♯ , ỹ♯ , z♯ ) provides a potential deployment plan which is feasible to the original design constraints and to a relaxed version of the operation constraints. However, x ♯ may be infeasible with respect to some of the original operation constraints containing integer operation variables. In order to discuss this point in detail, we need to come back to the initial formulation of CP.

Indeed, most of the feasibility issues encountered when solving the operation sub-problems for a given potential deployment plan found at the upper level come from the impossibility to comply with a specific set of operation constraints: the minimum production constraints. Namely, when the integer operation variable S t,m representing the number of active devices of type m at time step t is relaxed to a continuous variable, the corresponding continuous energy production variable P out t,n,m may take any value in the range [0, P max n,m ] rather than in the range {0} ∪ [P min n,m , P max n,m ]. The same difficulty arises when the working range of a set of devices is linked to the decision variable MODE t,n,m representing the number of devices operating under a given mode n. In other words, a solution of SRP complies with the original maximum production constraints but does not have to respect the minimum production constraints. As a result, the potential deployment plans found when solving SRP at the upper level often involve only high capacity devices which are efficient at producing a large amount of supply commodity but are not able to meet a small but non-zero demand due to the minimum power restrictions.

In this subsection, we thus investigate a set of valid inequalities to tighten the formulation of SRP and reduce the gap between CP and its relaxation SRP.

These valid inequalities focus on the supply commodities that cannot be sold to outside buyers, i.e. to the commodities in C S \ C S,Sell . The demand for supply commodity c = (n, s) can be satisfied by conversion and storage devices. Let us denote by c ′ = (n ′ , s ′ ) the input commodity of these storage devices. For example, in the trigeneration system in Figure 3.4, the supply commodity c = (HEAT, 2) has an associated intermediate stored commodity c ′ = (HEAT, 1) and the supply commodity c = (COLD, 1) has an associated intermediate stored commodity c ′ = (ICE, 1). We define the set of conversion device types which can produce either supply commodity c or its corresponding intermediate stored commodity

c ′ as M(c) = {(p, l) ∈ M CONV |p ∈ P in f low c ∪ P in f low c ′ } .
We introduce the concept of oversized devices. A device is defined as oversized with respect to a supply commodity c = (n, s) for a representative day (ϕ, d) if it will not be able to single-handedly meet the demand for this commodity for all time steps h ∈ {0, ..., H -1} of this day due to the restriction on its minimum power production. Namely, an oversized device cannot be used to satisfy the demand for c either because its minimum hourly production per hour of c is higher than the demand value of one hour in the day, or because its minimum production per hour of c ′ is higher than the total daily demand. We thus define two sets of oversized devices:

• the subset of M(c) of device types which do not produce c or whose minimum production per hour of c is too large to allow them to produce the lowest hourly demand for c = (n, s) during day (ϕ, d), i.e. P min m,n > min h=0,..,H-1 Dem ϕ,d,h,n,s . This subset is denoted as M O,1 ϕ,d,c . • the subset of M(c) of devices types which do not produce c ′ or whose minimum production per hour of c ′ is too large to allow them to produce the total daily demand of day (ϕ, d), i.e. P min m,n ′ > ∑ H-1 h=0 Dem ϕ,d,h,n,s . This subset is denoted as M O,2 ϕ,d,c . We define the set of oversized devices as the intersection of the oversized devices of these two types

M O ϕ,d,c = M O,1 ϕ,d,c ∩ M O,2 ϕ,d,c
. Clearly, a system layout at phase ϕ comprising devices which are all in the set M O ϕ,d,c will not be able to satisfy the demand for the supply commodity c for day d ∈ D ϕ while simultaneously complying with all the minimum power production constraints. Therefore, any feasible system layout at phase ϕ must involve at least one non oversized device, i.e. one device in M(c) not belonging to M O ϕ,d,c . Moreover, we note that, when expressing theses constraints, there is no need to consider all days in D ϕ . Namely, thanks to the way extreme days are selected (see Subsection 4.2.2), it is enough to focus, for each phase, on the extreme day displaying the smallest hourly demand and on the one displaying the smallest total daily demand. The constraints relative to the other days will be redundant with the ones relative to these two extreme days.

Let D low ϕ ⊂ D ϕ be the set of these two extreme days. We add the following valid inequalities to the formulation of SRP in a pre-treatment step, i.e. before starting the upper level branch-and-cut search.

∑ m∈M(c)\M

O ϕ,d ϕ ∑ φ=1 SD φ,m ≥ 1 ∀d ∈ D low ϕ , ∀ϕ (5.29)
By adding Constraint (5.29), we ensure that any system layout found in the upper level problem will contain at least one non-oversized device corresponding to the supply commodity c to deal with the energy demand with small values.

Moreover, due the definition of the set M O ϕ,d , we may encounter system layouts satisfying Constraints (5.29) but containing no devices in the set m ∈ M(c) \ M O,1 ϕ,d .

In this case, the only way the LMES may satisfy a small daily demand of c consists in producing c ′ during one (or several) hour(s) and storing it in an energy storage device. This implies that, if ∑ m∈M(c)\M O,1 ϕ,d ∑ ϕ φ=1 SD φ,m = 0, at least one unit of discrete energy storage capacity should be installed. This leads to a second set of valid inequalities:

ϕ ∑ φ=1 ∑ m ′ ∈M ′ (c ′ ) SD φ,m ′ ≥ 1 - ∑ m∈M(c)\M O,1 ϕ,d ϕ ∑ φ=1 SD φ,m ∀d ∈ D low ϕ , ∀ϕ (5.30)
The set M ′ (c ′ ) denote the set of storage devices m ′ = (p, l) such that p ∈ P STO ∩ P out f low c ′ , i.e. the set of storage devices able to store commododity c ′ . Constraints (5.30) state that if a system layout contains only oversized devices of Type 1, it should have at least one unit of discrete energy storage capacity corresponding to the stored commodity c ′ .

These valid inequalities will contribute in improving the efficiency of the hierarchical decomposition algorithm by a priori excluding some infeasible potential deployment plans, thus saving the computational effort needed to determine that they are infeasible.

Summary of the improved hierarchical decomposition algorithm

Before presenting our numerical results, we provide a summarized description of the proposed extended hierarchical decomposition algorithm: see Algorithm 2. Moreover, Figure 5.2 provides a graphical description of the improved hierarchical decomposition algorithm as applied to our multi-phase deployment problem.

Computational experiments

We now focus on assessing the performance of the hierarchical decomposition algorithm, both in its initial and extended variants, at solving the problem of optimally designing an LMES over a multi-phase investment plan. This assessment is carried out by using the set of instances based on our three case studies presented in Subsection 4.4.1.

Each instance is solved by the original hierarchical decomposition algorithm presented in [START_REF] Yokoyama | Optimization of energy supply systems by MILP branch and bound method in consideration of hierarchical relationship between design and operation[END_REF] and by the extended hierarchical decomposition discussed above. The experimental setup was the same as the one described at the beginning of Subsection 4.4.2. In both variants of the hierarchical decomposition algorithm, the upper level master problem was solved by CPLEX 12.8 solver. Incumbent Callbacks were used to control the flow of information between the upper and lower problems and Usercut Callbacks were used to set the upper cutoff value to the current value of Z best CP . Due to the presence of a Usercut Callback, which is classified as a "Control Callback", parallel MILP solving is turned off and CPLEX 12.8 solver uses a single thread to solve the upper level master problem. We imposed a time limit of 7200s.

Similar to what was done in Subsection 4.4.2, we provide, for each instance and each considered solution approach, the integrality gap Gap LP , the computation time and the residual gap Gap MIP . Moreover, we also report: Algorithm 2: Extended Hierarchical Decomposition Data: Parameters of commodities and technologies, construction phases and discount rate Result: The optimal system deployment plan and the minimum total cost 1 For each supply commodity in C S \ C • #IFS SRP : the number of integer feasible solutions of SRP explored during the upper level branch-and-cut search before the algorithm converges or the time limit is reached,

• #OP: the number of lower-level operation sub-problems solved over the course of the algorithm.

Results from Tables 5.1-5.3 first show that, in general, the initial hierarchical decomposition algorithm does not seem to be able to outperform the direct resolution by CPLEX 12.8. Namely, over the 19 considered instances, the average residual gap (Gap MIP ) is 1.59% when using this algorithm as compared to 0.40% when directly solving the whole problem as an MILP.

In contrast, the extended hierarchical decomposition algorithm provides a guaranteed optimal solution of the MILP problem formulated in Chapter 4 for all instances of City A project and for 4 out of 5 instances for City B project. Note that a near-optimal solution with a residual gap of 0.07% is obtained for the last instance related to City B project. Furthermore, the average computation time over the 12 instances related to City A or City B projects is significantly reduced from more than 7200s when using a direct resolution by CPLEX 12.8 or the initial hierarchical algorithm to 2949s when using the extended hierarchical decomposition algorithm.

This improved computational efficiency may be explained by three main reasons. First, the addition of valid inequalities (5.29)-(5.30) to SRP leads to a reduction of the integrality gap. Namely, over the 12 instances, the average value of Gap LP is reduced from 3.32% with the initial formulation of SRP to 2.56% with the formulation of SRP strengthened by valid inequalities (5.29)-(5.30). Even if this gap reduction is relatively small, it is obtained by adding a very limited number of constraints (4Φ) to the formulation and thus does not negatively impact the computation time needed to solve the linear relaxation of SRP at each node of the upper level B&C search tree. Moreover, adding these valid inequalities to the formulation enables us to a priori exclude from the upper level search space deployment plans which are feasible for SRP but not for CP and thus to save the computation time needed to determine that they are infeasible. Second, over the 12 considered instances, the average value of #IFS SRP , the number of integer feasible solutions of SRP explored during the upper level B&C search, is significantly reduced from 2467 with the initial hierarchical decomposition algorithm to 962 with the extended hierarchical decomposition algorithm. This shows the efficiency of the no-good cuts (5.28) at preventing the upper level B&C algorithm from uselessly considering multiple times a system layout which is known to be infeasible for a given investment phase. Finally, thanks to the combined used of the no-good cuts (5.28) and the recording of previously computed results in a hash table, the average number of lower-level operation sub-problems solved over the course of the hierarchical decomposition algorithm, #OP, is drastically reduced from 31709 with the initial variant of the algorithm to 894 with the extended variant. Even if each of these sub-problems is a small MILP which can be solved without any numerical difficulty, dividing by a factor of 35 the number of sub-problems solved over the course of the algorithm positively impacts the overall computation time.

However, the extended hierarchical decomposition algorithm seems to have a poor performance on the instances based on the trigeneration system located in City C. Namely, even if the integrality gap was significantly reduced by adding valid inequalities (5.29)-(5.30), out-of-memory issues were encountered while solving all the instances. Moreover, the average residual gap obtained with the extended hierarchical decomposition algorithm (Gap MIP = 0.53%) is larger than the one obtained with the direct resolution approach (Gap MIP = 0.37%). This difficulty may come from the large number of design variables. Namely, the total number of design variables, excluding the binary expansion variables, is 28 (resp. 25) for the DCS project of City A (resp. City B) , while there is a total of 36 design variables for the trigeneration project. This larger number of variables leads to a larger feasible space for SRP to be explored by the B&C algorithm and thus to a larger computational effort. A second reason might be the fact that the lower bound provided by Z SRP is of lower quality for this set of instances than for the City A and City B instances, which forces the B&C algorithm to explore more nodes before converging.

Conclusion

In this chapter, we first highlighted the bi-level structure of the MILP problem formulated in Chapter 4. This special structure allows to split the original large-scale MILP problem into a set of small-size MILP problems once the deployment plan of the system is determined.

We then applied the previously published hierarchical decomposition algorithm to the large-scale MILP problem formulated in Chapter 4. Note that the hierarchical decomposition algorithm requires that the upper-level variables are all discrete and that the original MILP model introduced in Chapter 4 meets this requirement. The hierarchical decomposition algorithm takes as upper-level master problem the semi-relaxed problem (SRP) which is obtained by relaxing all lower-level discrete variables of the original complete problem (CP). The resolution of the upper-level problem is carried out by the built-in B&C algorithm of a commercial mathematical solver. This B&C algorithm is however customized through the use of callbacks: at each node where an integer feasible solution of SRP, i.e. a potential deployment plan, is found, we solve a series of operation sub-problems (OP) to find out the real objective value corresponding to the current deployment plan. All integer solutions at the upper-level are rejected to guarantee the optimal convergence of the algorithm. We also proposed two ways to further improve the computational efficiency of the hierarchical decomposition method. The first improvement is related to the multi-phase nature of our design problem. In order to avoid repetitively solving several times the same sub-problems, we store in a hash table the optimal value of all feasible solved sub-problems and retrieve this value to reuse it when needed over the course of the algorithm. In case a sub-problem is found to be infeasible, we also add a no-good cut to exclude the current system layout for the corresponding investment phase from the feasible space of SRP. The second improvement consists in adding valid inequalities into the formulation of SRP. Namely, due to the relaxation of discrete lower-level variables, the SRP problem does not consider the minimum load rate of devices and tends to favor the installation of large-capacity devices. The proposed valid inequalities aim at preventing the SRP from choosing only oversized devices and thus contribute in reducing the number of infeasible deployment plans explored by the upper-level B&C algorithm.

|D

Both the initial and extended hierarchical decomposition algorithms are tested on the instances based on the three case studies introduced in Chapter 4. One one hand, the numerical results show that the initial hierarchical decomposition method is on average less efficient than the direct resolution of CP by CPLEX 12.8 solver. On the other hand, the extended hierarchical decomposition algorithm significantly outperforms the two other methods on the instances related to the two DCS projects, providing optimal (or near optimal) solutions for all instances within the time limit. However, out-of-memory issues were still encountered for the instances related to the trigeneration project, even the ones involving a small number of representative days.

This motivated us to explore another resolution strategy based on a generalized Benders' decomposition algorithm taking advantage of a specific property of the sub-problems of the MILP problem. This solution approach will be presented in Chapter 6.

the strong duality properties holding for linear and convex programs. In what follows, we discuss how this generalized Benders' decomposition framework may be applied to our LMES optimal design problem.

To this aim, we start by recalling the compact reformulation of the original problem, denoted by CP in Chapter 5.

Z * CP = min f 0 (x) + Φ ∑ ϕ=1 ∑ d∈D ϕ f ϕ,d (y ϕ,d , z ϕ,d ) (6.1) h(x) ≤ 0 (6.2) g ϕ,d (x ϕ , y ϕ,d , z ϕ,d ) ≤ 0 ∀(ϕ, d) ∈ D (6.3) l ϕ,d (y ϕ,d , z ϕ,d ) ≤ 0 ∀(ϕ, d) ∈ D (6.4)
x ∈ Z ν (6.5)

y ϕ,d ∈ R µ ∀(ϕ, d) ∈ D (6.6) z ϕ,d ∈ Z λ ∀(ϕ, d) ∈ D (6.7)
As mentioned in [START_REF] Bolusani | A Framework for Generalized Benders' Decomposition and Its Application to Multilevel Optimization[END_REF], a Benders' decomposition approach can be seen as a reformulation technique in which the original optimization problem is reformulated using only a subset of the original variables. This subset of variables are referred to as the first-stage variables whereas the variables not included in the reformulation are called the second-stage variables. Thus, a Benders' decomposition method essentially carries out a projection of the original optimization problem onto the subspace corresponding to the first-stage variables.

In order to explain how this reformulation/projection technique can be applied on CP, we first introduce some additional notation. Let X = {x ∈ Z ν |h(x) ≤ 0} denote the set of deployment plans feasible with respect to the design constraints (6.2) and X ϕ = proj x ϕ (X ) the projection of X onto the space of the design variables relative to phase ϕ. Let S ϕ = {(y ϕ , z ϕ ) ∈ R µ|D ϕ | × Z λ|D ϕ | |l ϕ,d (y ϕ,d , z ϕ,d ) ≤ 0, ∀d ∈ D ϕ } denote the set of operation schedules relative to days d ∈ D ϕ feasible with respect to the operation constraints (6.4). Finally, F ϕ = {(x ϕ , y ϕ , z ϕ ) ∈ X ϕ × S ϕ |g ϕ,d (x ϕ , y ϕ,d , z ϕ,d ) ≤ 0, , ∀d ∈ D ϕ } represents the set of system layouts and operation schedules relative to phase ϕ complying with the design, coupling and operation constraints (6.2)- (6.4).

Using this notation, CP can be equivalently formulated as:

Z * CP = min x∈X f 0 (x) + Φ ∑ ϕ=1 min (y ϕ ,z ϕ )∈S ϕ ∑ (ϕ,d)∈D ϕ f ϕ,d (y ϕ,d , z ϕ,d )|(x ϕ , y ϕ , z ϕ ) ∈ F ϕ (6.8)
We now replace each inner minimization problem by a function. This allows us to reformulate CP using only the first-stage design variable x:

Z * CP = min x∈X f 0 (x) + Φ ∑ ϕ=1 Z ϕ (x ϕ ) (6.9)
Function Z ϕ : R ν ϕ → R ∪ {+∞} is defined as:

Z ϕ (x ϕ ) =    min (y ϕ ,z ϕ )∈S ϕ ∑ d∈D ϕ f ϕ,d (y ϕ,d , z ϕ,d )|(x ϕ , y ϕ , z ϕ ) ∈ F ϕ if x ϕ ∈ proj x ϕ (F ϕ )
+∞ otherwise where proj x ϕ (F ϕ ) ⊂ X ϕ is the projection of F ϕ onto the space of the first-stage variables x ϕ . proj x ϕ (F ϕ ) can be intuitively seen as the set of system layouts for phase ϕ, x ϕ , such that there is at least one feasible operation schedule for each day d ∈ D ϕ . Z ϕ is thus a function that returns the sum of the cost of the best operation schedule that can be found for each day d ∈ D ϕ under the system layout x ϕ if one exists. In fact, Z ϕ (x ϕ ) is the sum, over all days d ∈ D ϕ , of the optimal value of the operation sub-problems OP ϕ,d (x ϕ ) introduced in Section 5.2. Moreover, Z ϕ prevents any first-stage design solution x ϕ that is not feasible with respect to one or several days d ∈ D ϕ , i.e. which is not in proj x ϕ (F ϕ ), from being considered.

In principle, solving Problem (6.9) would provide us with an optimal solution of CP. However, this would require a closed-form expression of each function Z ϕ . As noted in [START_REF] Bolusani | A Framework for Generalized Benders' Decomposition and Its Application to Multilevel Optimization[END_REF], even when such a closed form exists in theory, its description is typically of exponential size and computationally impractical.

The basic idea of a Benders' decomposition approach consists in overcoming this difficulty by replacing each function Z ϕ by a dual bounding function Z ϕ to obtain a relaxation of Problem (6.9) called the master problem.

A dual bounding function

Z ϕ : R ν ϕ → R ∪ {+∞} is defined as a function such that Z ϕ (x ϕ ) ≤ Z ϕ (x ϕ ) for all x ϕ ∈ R ν ϕ . This dual function is called strong at xϕ ∈ R ν ϕ if Z ϕ ( xϕ ) = Z ϕ ( xϕ ).
The master problem can thus be written as:

min x∈X f 0 (x) + Φ ∑ ϕ=1 Z ϕ (x ϕ ) (6.10)
Problem (6.10) is a relaxation of Problem (6.9) (and consequently a relaxation of the initial problem CP) and thus provides a lower bound of Z * CP . The aim of a generalized Benders' decomposition algorithm will be to iteratively improve the formulation of the master problem by progressively strengthening the dual bounding function Z ϕ used for each phase ϕ. Basically, at each iteration k, we solve the master problem using the current expression of the dual bounding functions, yielding a lower bound LB k of Z * CP and a design solution denoted by x k . We then evaluate Z ϕ (x k ϕ ), for each phase ϕ, by solving a series of operation sub-problems

OP ϕ,d (x k ϕ ), obtain an upper bound UB k = f 0 (x k ) + ∑ Φ ϕ=1 Z ϕ (x k ϕ ) of Z *
CP and use the information obtained during this process to update each dual bounding function Z ϕ so as to strengthen it. The algorithm terminates when LB k = UB k .

Hooker and Ottosson [START_REF] Hooker | Logic-based Benders decomposition[END_REF] showed that the finite convergence towards an optimal solution of such a generalized Benders' decomposition algorithm is guaranteed as long as two conditions are respected. The first one is that the approximation of the dual bounding function used at iteration k for each phase ϕ is strong at each of the previously considered iterates, i.e. is such that Z ϕ (x

j ϕ ) = Z ϕ (x j ϕ ) for all iterations j in 1...k. The second one is that each set proj x ϕ (F ϕ ) is finite, i.e. that |proj x ϕ (F ϕ )| < ∞.
The first condition can be more easily satisfied by defining Z ϕ as the maximum of the lower bounding functions obtained at each previous iteration, i.e. by defining Z ϕ at iteration k as:

Z ϕ (x ϕ ) = max 1≤j≤k Z j ϕ (x ϕ ) (6.11)
where Z j ϕ is the dual bounding function obtained at iteration j ≤ k. In this case, if Z j ϕ is strong at x j ϕ , the dual bounding function used at a latter iteration k ≥ j, Z ϕ , will also be strong at x j ϕ . The master problem (6.10) can be reformulated using auxiliary variables η ϕ to eliminate the maximum operator involved in (6.11). We thus obtain at iteration k the following master problem denoted by MP k .

min Z MP k = f 0 (x) + ϕ ∑ ϕ=1 η ϕ (6.12)
x ∈ X (6.13)

η ϕ ≥ Z j ϕ (x ϕ ) ∀ϕ ∈ {1, • • • , Φ}, ∀j ∈ {1, • • • , k} (6.14) F j ϕ (x ϕ ) ≤ 0 ∀ϕ ∈ {1, • • • , Φ}, ∀j ∈ {1, • • • , k} (6.15
)

η ϕ ∈ R ∀ϕ ∈ {1, • • • , Φ} (6.16)
In this formulation, Constraints (6.14) can be understood as Benders' optimality cuts. Moreover, in order to avoid handling infinite values of the lower bounding function Z ϕ , we may add a series of constraints to explicitly exclude first-stage variables x such that x ϕ / ∈ proj x ϕ (F ϕ,d ) for some phase ϕ. This may be done by adding some constraints called Benders' feasibility cuts. This is the purpose of Constraints (6.15) in which F j ϕ (x ϕ ) denotes a general mathematical expression involving elements of x ϕ .

Moreover, in our problem, all the first-stage decision variables are bounded integer variables: see the design constraints (4.1)-(4.2) introduced in Chapter 4. X ϕ and as a consequence proj x ϕ (F ϕ ) ⊂ X ϕ are therefore finite sets for all phases ϕ.

Hence, the main remaining challenge towards ensuring that the generalized Benders' decomposition will finitely converge towards an optimal solution is to build, for each phase ϕ, a strong lower bounding function at each iteration k, i.e. to build a function

Z k ϕ such that Z k ϕ (x ϕ ) ≤ Z ϕ (x ϕ ), ∀x ϕ ∈ R ν ϕ and Z k ϕ (x k ϕ ) = Z ϕ (x k ϕ )
. This is the purpose of the next section.

Strong dual bounding functions

We now investigate the development of strong dual bounding functions for our problem. This development exploits the specific structure of the coupling constraints (6.3) in CP.

We start by defining a partial order between the integer vectors describing the system layout at phase ϕ. Thus, a system layout described by vector x ♯ ϕ is defined to be smaller than the system layout described by vector x

• ϕ,d , Z ♯ ϕ,d ≥ Z • ϕ,d . Moreover, if OP ϕ,d (x • ϕ ) is infeasible, OP ϕ,d (x ♯
ϕ ) is also infeasible. Note that this mathematical reasoning can be intuitively related to the fact that a system layout x ♯ ϕ which has more components of each type (i.e. which has more energy conversion devices of each type, a larger storage capacity, a larger contracted power....) than the system layout x • ϕ can be operated using the same operation schedules as the ones used with the system layout x • ϕ (and maybe with less expensive operation schedules). This implies that the operation cost of x ♯ ϕ cannot be larger than the one of x • ϕ . This property can be illustrated by an example based on the case study related to the DCS project located in City B. We consider the operation cost for Day 143 of Phase 4 of this project. The number of single-mode chillers of type (SMEC, 1) is fixed to 0, the ice storage capacity to 171MWh and the contracted maximum power to 15MW. Figure 6.1 shows how the operation cost for this day, Z 4,143 , varies as a function of the design vector x 4 = (x 4,1 , x 4,2 ) where x 4,1 is the number of chillers of type (SMEC, 2) and x 4,2 the number of chillers of type (DMEC, 1) present in the DCS at the beginning of Phase 4. The z-axis indicates the value of the operation cost Z 4,143 (x 4 ): a value of Z 4,143 (x 4 ) above 550000 corresponds to an 'infinite value' indicating that the system layout x 4 is infeasible for day (4,143). Note that, in order to obtain a more readable figure, x 4 was considered as a continuous vector when drawing the graphical representation of Z 4,143 but the scheduling variables z 4,143 were kept integer when solving OP 4,143 (x 4 ).

We first observe that OP 4,143 (x 4 ) is infeasible for small values of x 4 (see the region colored in yellow on the figure). Moreover, if OP 4,143 (x ♯ 4 ) is infeasible for a given design vector x ♯ 4 , it is infeasible for all smaller design vectors: see e.g. how vector (1,3) and all smaller vectors, i.e. vectors (0,3) (0,2), (0,1), (1, 2), (1,1), (1,0) and (0,0), are infeasible. We also note that Z 4,143 is a piece-wise constant function and that for a given value of x 4,1 (resp. x 4,2 ), Z 4,143 is non-increasing with respect to x 4,2 (resp. x 4,1 ). Finally, we observe that if OP

∈ D ϕ , Z ϕ,d (x ♯ ϕ ) ≤ Z ϕ,d (x ϕ ) for any x ϕ such that x ϕ ⪯ x ♯ ϕ .
In other words, Z ϕ,d (x ♯ ϕ ) is a lower bound on the operation cost of day d ∈ D ϕ as long as x ϕ is smaller than x ♯ ϕ . This allows us to define a lower bounding function for the total operation cost at phase ϕ, Z ϕ , as follows.

Z ϕ (x ϕ ) = Z ϕ (x ♯ ϕ )1 x ϕ ⪯x ♯ ϕ (6.17)
where 1 x ϕ ⪯x k ϕ is the indicator function defined by:

1 x ϕ ⪯x ♯ ϕ = 1 if x ϕ ⪯ x ♯ ϕ 0 otherwise (6.18) We have Z ϕ (x ϕ ) = Z ϕ (x ♯ ϕ ) ≤ Z ϕ (x ϕ ) for any x ϕ smaller that x ♯ ϕ and Z ϕ (x ϕ ) = 0 ≤ Z ϕ (x ϕ ) for any x ϕ greater or incomparable to x ♯ ϕ . Thus Z ϕ (x ϕ ) ≤ Z ϕ (x ϕ ) for all x ϕ ∈ R ν ϕ . Z ϕ is thus a lower bounding function of Z ϕ . Moreover, we have Z ϕ (x ♯ ϕ ) = Z ϕ (x ♯ ϕ ) so that Z ϕ is strong at x ♯ ϕ . Furthermore, let x ♯ ϕ /
∈ proj x ϕ (F ϕ ) be a system layout infeasible for some day d ∈ D ϕ . As discussed above, we know than any x ϕ such that x ϕ is smaller that x ♯ ϕ will also be an infeasible system layout for this day. We can thus exclude any x ϕ smaller that x ♯ ϕ from the feasible space of the master problem by adding the following equality:

1 x ϕ ⪯x ♯ ϕ = 0 (6.19)
Thus, the master problem MP k to be solved at iteration k of a generalized Benders' decomposition algorithm may be rewritten as:

min Z MP k = f 0 (x) + ϕ ∑ ϕ=1 η ϕ (6.20)
x ∈ X (6.21)

η ϕ ≥ Z ϕ (x j ϕ )1 x ϕ ⪯x j ϕ ∀ϕ ∈ {1, • • • , Φ}, ∀j ∈ {1, • • • , k|x j ϕ ∈ proj x ϕ (F ϕ )} (6.22) 1 x ϕ ⪯x j ϕ = 0, ∀ϕ ∈ {1, • • • , Φ}, ∀j ∈ {1, • • • , k|x j ϕ / ∈ proj x ϕ (F ϕ )} (6.23) η ϕ ∈ R ∀ϕ ∈ {1, • • • , Φ} (6.24)
In this formulation, Constraints (6.22) are optimality cuts generated at each iteration j to approximate the total operation cost at phase ϕ when the iterate x j ϕ belongs to proj x ϕ (F ϕ ). In case x j ϕ / ∈ proj x ϕ (F ϕ ), i.e. in case there is no feasible operation schedule for some day d ∈ D ϕ with the system layout x j ϕ , we add a feasibility cut (6.23) to forbid x j ϕ .

MILP reformulation of the master problem

Problem (6.20)-(6.24), as such, cannot be directly solved by an MILP solver due to the presence of a set of non-linear indicator functions 1 x ϕ ⪯x j ϕ , j = 1, • • • , k. We thus investigate in what follows a MILP reformulation of the master problem MP k .

To this aim, we introduce two sets of binary variables. First, for each j = 1, • • • , k and each phase ϕ, the binary variable β j ϕ is defined as

β j ϕ = 1 if x ϕ ⪯ x j ϕ and β j ϕ = 0 otherwise. We thus have β j ϕ = 1 x ϕ ⪯x j ϕ . Moreover, for each j = 1, • • • , k, each phase ϕ and each element i = 1, • • • , ν ϕ of x j ϕ , the binary variable α j ϕ,i is defined by α j ϕ,i = 1 if x ϕ,i > x j ϕ,i and α j ϕ,i = 0 otherwise. We thus have β j ϕ = 1 iif α j ϕ,i = 0 for all i = 1, • • • , ν ϕ .
Thanks to these binary variables, the master problem MP k can be reformulated as follows:

min Z MP k = f 0 (x) + ϕ ∑ ϕ=1 η ϕ (6.25)
x ∈ X (6.26)

η ϕ ≥ Z ϕ (x j ϕ )β j ϕ ∀ϕ, ∀j ∈ {1, • • • , k|x j ϕ ∈ proj x ϕ (F ϕ,d )} (6.27) β j ϕ = 0 ∀ϕ, ∀j ∈ {1, • • • , k|x j ϕ / ∈ proj x ϕ (F ϕ,d )} (6.28) α j ϕ,i ≤ x ϕ,i x j ϕ,i + 1 ∀ϕ, ∀j ∈ {1, • • • , k}, ∀i ∈ {1, • • • , ν ϕ } (6.29) α j ϕ,i ≥ x ϕ,i -x j ϕ,i x ϕ,i -x j ϕ,i ∀ϕ, ∀j ∈ {1, • • • , k}, ∀i ∈ {1, • • • , ν ϕ } (6.30) β j ϕ ≥ 1 - ν ϕ ∑ i=1 α j ϕ,i ∀ϕ, ∀j ∈ {1, • • • , k} (6.31) β j ϕ ≤ ∑ ν ϕ i=1 (1 -α j ϕ,i ) ν ϕ ∀ϕ, ∀j ∈ {1, • • • , k} (6.32) η ϕ ∈ R ∀ϕ (6.33) β j ϕ ∈ {0, 1} ∀ϕ, ∀j ∈ {1, • • • , k} (6.34)
Constraints (6.27)-(6.28) are the optimality and feasibility cuts (6.22)-(6.23) in which the indicator function 1 x ϕ ⪯x j ϕ has been replaced by the binary variable β j ϕ . Constraints (6.29)-(6.30) link the value of α j ϕ,i to the one of element i of vector x j ϕ . They ensure that α j ϕ,i = 0 if x ϕ,i ≤ x j ϕ,i (see Constraints (6.29)) and α j ϕ,i = 1 if x ϕ,i > x j ϕ,i (see Constraints (6.30) in which x ϕ,i denotes an upper bound on x ϕ,i ). Finally, Constraints (6.31)-(6.32) link the value of β j ϕ to the value of variables α j ϕ,i . Thus, Constraint (6.31) forces β j ϕ to be equal to 1 when all variables α j ϕ,i , ∀i ∈ {1, • • • , ν ϕ }, are equal to 0, i.e. when all elements of x ϕ are smaller that their counterpart in x j ϕ . Similarly, Constraint (6.32) ensures that β j ϕ is equal to 0 as soon as there is one element of x ϕ which is greater that that its counterpart in x j ϕ , i.e. as soon as one of the variables α j ϕ,i , i ∈ {1, • • • , ν ϕ } is equal to 0. It should be mentioned that the upper bound x ϕ,i should numerically be set to a value large enough to ensure that the denominator of the right-hand-side of Constraints (6.30) remains strictly positive and to avoid the division-by-zero issue.

Generalized Benders' decomposition algorithm

Sections 6.2 to 6.4 described a framework for developing a generalized Benders' decomposition approach applied to our problem and provided a set of lower bounding functions displaying the sufficient mathematical properties to ensure the finite and optimal convergence of this algorithm. Overall, the outline of the algorithm can be described by Algorithm 3. However, a direct implementation of Algorithm 3 showed some numerical inefficiencies.

if OP ϕ,d (x k ϕ ) is feasible then 14 set Z ϕ (x k ϕ ) to Z ϕ (x k ϕ ) + Z ϕ,
Set UB k = min{ f 0 (x k ) + ∑ Φ ϕ=1 Z ϕ (x k ϕ ), UB k-1 }; 29 if |UB k -LB k | < ϵ
First, the solution of MP k found at Line 4 of Algorithm 3 may be far from being feasible (never mind optimal) for the operation sub-problems, especially in the first iterations of the generalized Benders' decomposition algorithm. In particular, we noted that x k is sometimes neither feasible nor optimal with respect to the linear relaxation of the operation sub-problems denoted by OP ϕ,d (x k ϕ ). As a consequence, when running Algorithm 3, we sometimes solve a series of operation sub-problems OP ϕ,d (x k ϕ ) as MILPs only to find out that x k ϕ is infeasible with respect to the linear relaxation of one of this sub-problems. This computational effort could be avoided by first strengthening the formulation of MP k through the use of 'classical' Benders' cuts generated using the strong duality property holding for the linear programs OP ϕ,d (x k ϕ ). These cuts may be generated through a significantly smaller effort than the one required to generate cuts of type (6.27)-(6.28) as they only necessitate the resolution of a set of LPs rather than MILPs. This is why we propose to replace the resolution of MP k with an MILP solver by the resolution of another relaxation of CP, denoted by SRP k , through a classical Benders' decomposition algorithm implemented in an off-the-shelf MILP solver. Problem SRP k corresponds to Problem SRP introduced in Chapter 5, i.e. to a relaxation of CP in which all integer variables y are relaxed to continuous variables, in which the cuts of type (6.27)-(6.28) generated up to iteration k by the generalized Benders' decomposition algorithm have been added. This allows us to obtain, at each iteration of the generalized Benders' decomposition algorithm, first-stage solutions x k , k = 1...K, of better quality and to solve operation sub-problems OP ϕ,d (x k ϕ ) as MILPs, which requires a significant computational effort, only to evaluate potentially 'interesting' deployment plans.

Second, similar to what was observed for the hierarchical decomposition algorithm, we often repetitively solve the same lower level sub-problems over the course of the algorithm. This repetition is due to the multi-phase nature of our problem. Namely, even if two deployment plans found by solving the master problem differ from each other, they may contain the same system layout for a subset of construction phases. This is why we record in a hash table the optimal value Z ϕ (x k ϕ ) of the total operation cost relative to a given phase ϕ obtained with the system layout x k ϕ considered at iteration k and retrieve this cost at a later iteration j > k of the algorithm in which the same system layout x j ϕ = x k ϕ is used for phase ϕ. This leads to an extended variant of the generalized Benders' decomposition algorith described in Algorithm 4.

Note that solving SRP k using the classical Benders' decomposition algorithm embedded e.g. in CPLEX solver at the upper level amounts to adding to MP k a set of optimality and feasibility cuts based on the resolution of the dual of each relaxed operation sub-problem OP ϕ,d (•). We thus build a 'temporary' dual bounding function for each function Z ϕ , ϕ ∈ {1, • • • , Φ}, involved in Problem (6.9), which is most often not strong at iterate x k . However, these dual bounding functions are only used to improve the quality of iterate x k . Over the course of each iteration, we build a dual bounding function for each function Z ϕ which is strong at iterate x k . This ensures that Algorithm 4 has a finite and optimal convergence.

if OP ϕ,d (x k ϕ ) is feasible then 18 set Z ϕ (x k ϕ ) to Z ϕ (x k ϕ ) + Z ϕ,
Set UB k = min{ f 0 (x k ) + ∑ Φ ϕ=1 Z ϕ (x k ϕ ), UB k-1 }; 32 if |UB k -LB k | < ϵ

Computational experiments

We report in this section the results of the computational experiments carried out to assess the numerical performance of the extended variant of the generalized Benders' decomposition algorithm introduced in Section 6.5 and to compare it with the one of the direct resolution by an MILP solver and of the hierarchical decomposition algorithms investigated in Chapter 5. We then discuss the outcome of a postoptimization simulation study carried out to evaluate the quality of the deployment plans obtained for the two DCS projects.

Numerical efficiency of the extended generalized Benders' decomposition algorithm

The computational experiments were based on instances related to the three case studies introduced in Chapter 4. Each instance was solved with the extended variant of the generalized Benders' decomposition algorithm. The experimental setup was the same as the one described at the beginning of Subsection 4.4.2. We use the classical Benders' decomposition algorithm embedded in CPLEX 12.8 solver to solve SRP k at each iteration of the algorithm.

The results are shown in Tables 6.1 -6.4. We report, for each instance, the value of the residual gap Gap MIP , the computation time and the number of iterations #Iter carried out by the generalized Benders' algorithm.

We first observe from the results displayed in Tables 6.1, 6.2 and 6.4 that the modified generalized Benders' decomposition algorithm significantly outperforms both the direct resolution approach and the extended hierarchical decomposition for all considered instances. Namely, it was able to solve to optimality the 19 considered instances within the imposed time and memory limits whereas the extended hierarchical decomposition algorithm could only solve to optimality 11 out of these 19 instances and encountered out-of-memory issues for the 7 instances related to the trigeneration project. Moreover, the modified generalized Benders' decomposition algorithm has a significantly higher efficiency than the extended hierarchical decomposition algorithm. Thus, over the 12 instances corresponding to the DCS projects in Cities A and B, the average computational time is decreased from 2279s with the extended hierarchical decomposition to 628s with the generalized Benders' decomposition, which represents a reduction of the computation time by a factor larger than 3.6.

Furthermore, since the average computational time for the 5 instances related to the DCS project in City B was less than 1 minute (see Table 6.2), we created 6 additional larger instances for this project by increasing the number of typical days taken into account to represent the demand evolution over the year. Note that the instance with |D ϕ | = 165 takes into account the 165 days with a non-zero cooling demand in the initial time series and assigns a weight equal to 1 to each of them, it is thus equivalent to the initial problem considering the cooling demand of all 365 days of a year. We note that, although the computational time increases with the size of the instance, the generalized Benders' decomposition algorithm returns the optimal solution for these 6 new instances within the time limit.

The second observation deals with the number of iterations #Iter carried out by the generalized Benders' decomposition algorithm. It seems that, for a given case study, this number does not change when the instance size increases, i.e. when the number of typical days |D ϕ | increases. This means that the increase of the computation time observed when |D ϕ | increases mainly comes from the additional computational effort needed to solve the larger number of operation sub-problems. However, we note that the number of iterations #Iter varies significantly with the case study, e.g. from an average value of 4 iterations for City B project to an average value of 20 iterations for City C project. These variations may be first explained by the number of design variables. This one is namely smaller in City B projects (25 design variables) than in City C project (36 design variables). This means that the design feasible space is larger for City C project than for City B project, which may make it harder for the generalized Benders' decomposition algorithm to converge towards the optimal deployment plan. A second explanation might be related to the 'quality' of the deployment plans obtained while solving Problem SRP k with a classical Benders' decomposition algorithm. This quality depends heavily on the strength of formulation SRP k : if this formulation is strong, it is more likely to provide deployment plans which will be feasible (or even close to optimal) for the complete problem CP so that the generalized Benders' algorithm will have to carry out less iterations. This point should however be further investigated.

Discussion on the quality of the obtained deployment plans

Solving the MILP problem formulated in Chapter 4 to optimality does not guarantee that the obtained deployment plan is optimal with respect to the initial optimization problem. Indeed, in the MILP problem formulated in Chapter the consumption of all devices displaying non-linear performance curves is not computed exactly due to the use of a piece-wise linear approximation of these curves. A too rough estimate of the actual operation cost in the MILP problem may lead to poor design decisions and thus negatively impact the quality and relevancy of the proposed deployment plan. In this subsection, we thus seek to assess the quality of the deployment plans obtained with the proposed modeling and solution approach. We first discuss the impact of the number of selected days |D ϕ | on the deployment plans. Tables 6.5-6.8 describe respectively the best feasible deployment plan found for City A, City B and City C projects obtained by solving the MILP problem formulated in Chapter 4 with the generalized Benders' decomposition algorithm as a function of the number of representative days used to estimate the operation cost. We first provide the value of the objective function (4.35), i.e. the value of the total net present cost of the deployment plan computed by the MILP model, followed by the value of the actualized design cost and the value of the actualized operation cost. We then give, for each phase ϕ in {1, ..., Φ}, the number SD ϕ,m of devices of each type to be installed in each phase, the energy storage capacity Size step m SD ϕ,m to be built and the maximum instantaneous power C step (ELEC,1) Cdis ϕ,(ELEC,1) to be contracted. Results from Table 6.5 show that, for City A project, even if the total net present cost of the deployment plan slightly varies with the number of representative days, its actualized design cost and the design decisions to be made at each phase are the same in all the solutions obtained for |D ϕ | greater than 14 days. Similarly, we observe in Tables 6.6 -6.7 that, for City B project, the obtained deployment plans are the same for all values of |D ϕ |. Results from Table 6.8 also show that, for the trigeneration project, the optimal deployment plans are the same for |D ϕ | greater 6.9 for City A project and in Table 6.10 for City B project. We provide for each phase ϕ:

• Z optim ϕ : the value of the approximate operation cost in the objective function (4.35) of the optimization model (computed using a limited set of representative days and a piece-wise linear approximation of the chillers performance curves),

• Z sim ϕ : the value of the simulated operation cost obtained at the end of Step 2 of the simulation procedure (computing using the initial time series for the cooling demand and a piece-wise linear approximation of the chillers performance curves),

• Z actual ϕ : the value of the actual operation cost obtained at the end of Step 3 of the simulation procedure (computing using the initial time series for the cooling demand and the original non-linear chillers performance curves provided by the manufacturer),

• ϵ optim ϕ : the relative percentage difference between the value of the approximate operation cost computed by the optimization model and the value of the actual operation cost computed by the simulation procedure,

• ϵ sim ϕ : the relative percentage difference between the value of the simulated operation cost computed at Step 2 of the simulation procedure and the value of the actual operation cost computed by the simulation procedure.

The last line of each table provides the corresponding net present operation cost over the whole lifetime of the DCS.

We note from the results presented in Tables 6.9-6.10 that the approximation of the operation cost used in the optimization model is of very good quality. Namely, the relative percentage difference between Z optim ϕ and Z actual ϕ , ϵ optim ϕ , stays below 1% in all but one cases. The only large percentage difference is observed for Phase 1 of City B project. It is mainly explained by the fact that the demand in year 1 of this project is weak and that, in this case, it may be profitable to produce in one time-step an amount of ice covering the total cooling demand over the next two forthcoming days. In the optimization model, such a production strategy in which ice is kept in inventory for more than one day is not exploited due to the use of separate representative days. However, this has a limited impact on the relative percentage difference between the total approximated cost over all phases, Z optim = ∑ Φ ϕ=1 β ϕ Z optim ϕ

, and the total actual cost over all phases, Z actual ϕ = ∑ Φ ϕ=1 β ϕ Z actual ϕ , as the relative difference stays below 1% for both projects. We thus conclude that the use of a large number of representative days and of a 3-segment piece-wise linear approximation of the chillers performance curves leads to an estimation of the actual electric consumption of the system in the optimization model accurate enough to provide good-quality deployment plans. It would clearly be interesting to carry out such a post-optimization simulation study for the trigeneration project.

Conclusion

This chapter was devoted to the presentation of a generalized Benders' decomposition algorithm for the LMES design problem. We first presented in Section 6.2 the overall framework used to develop this algorithm. This framework basically consists in reformulating the problem using only design variables and in using a set of functions to represent the impact on the operation cost at a given investment phase of the design decisions. As there is no closed-loop expression for these functions, the algorithm will iteratively build an approximation, called a dual bounding function, for each of them.

The finite and optimal convergence of a generalized Benders' decomposition algorithm requires that these dual bounding functions are strong at all the iterates. We thus showed in Section 6.3 how such dual bounding functions may be used in our problem. This was achieved by exploiting the special structure of the constraints coupling design and operation variables in the operation sub-problems. We then explained in Section 6.4 how the resulting non-linear master problem may be reformulated as an MILP thanks to the introduction of a set of additional binary variables. Finally, Section 5.3 summarized the proposed generalized Benders' decomposition algorithm and highlighted the extension of the basic variant of the algorithm we introduced in order to improve its numerical efficiency.

Finally, Section 6.6 reported the results obtained on instances based on our three case studies. These results show that the generalized Benders' decomposition algorithm clearly outperforms the solution approaches investigated in the previous chapters. In particular, using this algorithm, we were able to solve to optimality all considered instances within the time and memory limits. Furthermore, we presented the results of a post-optimization simulation study carried out to estimate the quality of the deployment plans obtained with the proposed approach. These results show that, even if some approximations and simplifications were made when modeling the actual optimization problem, the obtained deployment plans are of good quality and may be recommended for a real-life implementation.

Chapter 7 Conclusion

Conclusion

This PhD thesis aimed to develop a numerical decision-aid tool based on mathematical optimization to help our industrial partner, EDF Chine, at optimally designing a Local Multi-Energy System (LMES). The LMES optimal design problem consists in choosing a combination of devices belonging to various energy conversion and storage technologies to compose a district energy system which can satisfy the forecasted energy demands of the local residents and minimize the total cost comprising the installation cost and the operation cost over the system's lifetime.

We provided in Chapter 2 an overview of the current state of the art on this problem. In particular, an analysis of the papers proposing mathematical programming models for the LMES design problem was carried out while focusing on four important features: the modeling of the design decisions (selection and sizing of the technologies), the representation of the performance curves of the conversion devices, the incorporation of an imposed minimum load for conversion devices and the presence of storage devices. This analysis enabled us to see that none of the previously published works could meet the needs of our industrial partner in terms of problem modeling. Namely, none of the reviewed works consider that building an LMES is a multi-phase project in which the system is gradually built and extended over the years. Moreover, the modeling of the design decisions is often oversimplified by assuming that a technology may be sized by choosing its capacity within a continuous range. This is not practically relevant as the sizing of a technology depends in fact of the models available in the catalog of the equipment manufacturer.

Chapter 3 was devoted to a detailed description of the LMES design problem. We started by introducing the two main components of an LMES: the commodities and the energy conversion and storage technologies. We explained how the LMES design problem uses as a starting point a description of the LMES superstructure which can be seen as a pre-selection of the technologies and corresponding devices considered for installation in the system. We then described a number of additional aspects to be considered when designing an LMES, in particular the fact that building an LMES requires a multi-phase deployment plan and that detailed operation schedules need to be built to accurately estimate the operation cost of the system.

In Chapter 4, we discussed the problem modeling and its formulation as an MILP. The proposed model extends the current state of the arts by simultaneously considering: (i) multiple construction phases, (ii) a realistic choice of the energy conversion and storage devices among the models available at the manufacturer, (iii) a careful selection of a number of representative days sufficiently large to obtain a good estimation of the actual operation cost of the system, (iv) a realistic representation of the actual functioning of the conversion devices in terms of minimum working load and partial-load efficiency. As taking into account all these aspects of the problem leads to the formulation of a huge mixed-integer non-linear program which, as such, is computationally intractable, we proposed a modeling approach aiming at reducing to some extent the size of the obtained mathematical program and at removing the non-linearities. This was achieved among others by using investment phases spanning several years, using a medium number of representative days to estimate the operation cost, building piece-wise linear approximations of the performance curves and exploiting the convexity of these curves. The proposed MILP formulation was tested on instances based on three case studies corresponding to two DCSs and a trigeneration system currently under development by EDF in China. However, significant numerical difficulties were encountered while trying to directly solve the proposed MILP formulation with CPLEX 12.8 solver.

We thus investigated the development of two decomposition-based approaches. These approaches exploit the bi-level nature of our problem, i.e. the fact that there is a natural hierarchy between the design decision variables used to determine the system deployment plan and the operation decision variables used to build operation schedules and to estimate the future operation cost of the system. This hierarchy translates into a special structure of the mathematical problem. Namely, once the design decisions determining the system layout at each investment phase are fixed, the problem decomposes into a set of independent operation sub-problems. Each of these sub-problems is a small-size MILP aiming at building an optimal operation schedule for a given representative day under the system layout determined for the corresponding phase.

The first algorithm was discussed in Chapter 5. It consists in an extension of the hierarchical decomposition algorithm presented by [START_REF] Yokoyama | Optimization of energy supply systems by MILP branch and bound method in consideration of hierarchical relationship between design and operation[END_REF]. This algorithm is based on the resolution of a semi-relaxed problem (SRP), i.e. of a relaxation of the original MILP in which all the integrality restrictions on the discrete operation variables are relaxed, by the built-in branch-and-cut algorithm of an MILP solver. Each time an integer feasible solution of SRP, i.e. a potential system deployment plan involving only integer design variables, is found in the B&C search tree, a sequence of operation sub-problems (OP) in which the operation variables are kept discrete when relevant is solved. This enables to first check the feasibility of the potential system deployment plan with respect to the operation constraints and second to compute the actual operation cost of this deployment plan. The nodes in the upper level B&C search tree are only pruned when their bound exceeds the value of the actual design and operation cost of the currently best known feasible deployment plan, so that the optimality of the final solution is guaranteed. We proposed a new extension of this hierarchical decomposition algorithm to make it more efficient at solving our multi-phase design problem. This was achieved in particular through the addition of single-phase no-good cuts into the upper-level problem to forbid deployment plans using a system layout found to be infeasible for a certain investment phase, through the recording of the solution value of previously solved operation sub-problems to avoid repetitive calculations and through the use of a small set of valid inequalities. Our numerical results showed that the proposed extended hierarchical decomposition algorithm significantly outperforms both the initial hierarchical decomposition algorithm and the direct resolution by a mathematical solver on the instances based on the two DCS projects. However, for the instances based on the trigeneration system, which has a more complex structure than the two DCSs, the extended hierarchical decomposition algorithm could not find the optimal solution within the computation time and memory limits. This is why we investigated the development of a generalized Benders' decomposition algorithm by using the generic framework provided in [START_REF] Bolusani | A Framework for Generalized Benders' Decomposition and Its Application to Multilevel Optimization[END_REF]. This algorithm can be seen as an extension of the conventional Benders' decomposition algorithm to the case where the second-stage sub-problems involve integer variables. It is based on the use of a master problem involving only design variables and of a set of dual bounding functions providing an estimation of the operation cost at each investment phase as a function of the system layout fixed for this phase by the design variables. The optimal and finite convergence of this algorithm requires among others that these dual bounding functions are strong. This was achieved in our case by proposing new optimality and feasibility cuts exploiting the special structure of the constraints linking the design and operation variables in our problem. Moreover, we developed an extended variant of this generalized Benders' decomposition algorithm in order to further improve its efficiency. This extended variant relies among others on the resolution of the semi-relaxed problem mentioned above by the classical Benders' decomposition algorithm embedded in an MILP solver. This allows us to improve the convergence speed of the algorithm by improving the quality of the potential deployment plan evaluated at each iteration of the generalized Benders' decomposition algorithm. Our computational results showed that the extended variant of the generalized Benders' decomposition algorithm clearly outperforms all the solution approaches investigated in the previous chapters. In particular, using this algorithm, we were able to solve to optimality all the considered instances, even the ones related to the trigeneration project, within the time and memory limits. However, solving the MILP problem to optimality does not guarantee that the obtained deployment plan is optimal with respect to the initial optimization problem as some approximations were made during the problem modeling. This is why we conducted a post-optimization simulation study to estimate the quality of the deployment plans obtained for the two DCS projects with the proposed approach. The outcome of this study was that, even if some approximations and simplifications were made when modeling the actual optimization problem, the obtained deployment plans are of good quality and may confidently be recommended for a real-life implementation.

Note that the work related to the hierarchical decomposition approach introduced in Chapter 5 was presented in several national conferences [Liu+19; Liu+20a; Liu+21b] and was published in the proceedings of an international conference [START_REF] Liu | A Hierarchical Decomposition Approach for the Optimal Design of a District Cooling System[END_REF]. A paper version was also submitted to an international journal [Liu+].

Research perspectives

The work presented in this manuscript raises several research perspectives that would be worth further investigation.

First, an interesting direction for further research might be to study the hybridisation of the hierarchical decomposition and generalized Benders' decomposition algorithms. For instance, the optimality and feasibility cuts of the generalized Benders' decomposition algorithm could be added during the resolution of the semirelaxed problem (SRP) by the upper level B&C algorithm of the hierarchical decomposition algorithm each time an integer feasible solution of SRP is found in the B&C search tree. This would avoid repetitively solving to optimality a sequence of variants of the semi-relaxed problem (SRP), differing from one another only by the addition of a small set of optimality and feasibility cuts, as is currently done in the generalized Benders' decomposition algorithm. However, this hybridisation is not straightforward and raises several difficulties. Namely, generating the optimality and feasibility cuts of the generalized Benders' decomposition algorithm requires the introduction of additional binary variables into the problem formulation during its resolution by a B&C algorithm (see Section 6.4). This means that the B&C algorithm would become a branch-and-price-and-cut algorithm. Furthermore, even if using the hybrid algorithm may reduce the time needed for obtaining integer feasible solutions of the original problem, it will involve however carrying out an accurate computation (through the resolution of a set of scheduling sub-problems) of the operation cost of potential deployment plans that may be of lesser quality than the optimal solution of the semi-relaxed sub-problem SRP. This may negatively impact the overall computation time of this hybrid approach.

Second, in terms of problem modeling, the MILP formulation we proposed for the LMES design problem relies on the assumption that the performance curves of all the conversion devices are linear or convex. This convexity has two advantages. First, it allows us to consider the epigraph of the curve rather than the curve itself so that the energy consumption of a device can be accurately computed as a function of its energy production without introducing auxiliary binary variables in the formulation. Second, the operation schedule can be built at the aggregate level for each set of identical devices (using a small number of integer variables) rather than at the detailed level for each individual device (using a larger number of binary variables). However, in case this assumption on the convexity of the performance curves does not hold, the formulation becomes more difficult to deal with. The piece-wise linear approximation method discussed in Subsection 4.2.3 can still be applied to approximate the non-convex performance curves. But the exact computation of the energy consumption of a device as a function of its energy production will require the use of formulation techniques based on auxiliary binary variables or Special Ordered Set (SOS) variables (see e.g. [START_REF] Sánchez-Martín | Gas network optimization: A comparison of piecewise linear models[END_REF]). Furthermore, it may not be possible anymore to build aggregate schedules for each set of identical devices. Namely, this aggregation is made possible by the property that all active identical devices should work at the same load when the performance curves are convex. This property is however not valid in the presence of non-convex performance curves. This would force us to build a detailed schedule for each individual device, which would require the introduction a large number of binary variables to represent the on/off status (and in some cases the operating mode) of each conversion device. Moreover, symmetry problems coming from the fact that there will exist many alternative optimal solutions that differ from one another only by the way each detailed operation schedule is assigned to one of the identical conversion devices may arise and lead to additional numerical difficulties.

Taking into account non-convex performance curves would thus complicate the resolution of the problem by the proposed decomposition algorithms, leading to a potential loss of numerical efficiency, but this resolution would still be possible. In contrast, some modeling extensions might make it impossible to use these decomposition algorithms. Namely, both the hierarchical decomposition and the generalized Benders' decomposition algorithms require that all design variables are discrete. In our case, this requirement could be rather easily met as the number of conversion devices to be installed should be chosen as an integer and the size of the energy storage to be built should be an integer multiple of a predefined small discrete capacity unit. However, in case a design variable must be considered as continuous, the optimal convergence of the two decomposition algorithms investigated in this work is not guaranteed anymore and we will have to solve the large-scale MILP problem directly with a mathematical solver. Furthermore, both decomposition algorithms are based on the key property that, once the design decisions are fixed, the problem decomposes into a set of independent scheduling sub-problems. This property might not hold anymore if the LMES to be designed involves long-term seasonal energy storage. This technology consists e.g. in producing and storing ice in winter to use it several months later to provide cold during summer. It thus differs from the shortterm daily energy storage considered in our work. Taking into account a seasonal energy storage means introducing in the MILP formulation an additional coupling (through a set of constraints) between the operation scheduling sub-problems. The development of a new solution approach might be required to handle this difficulty.

Finally, the current decision-aid tool uses as input data forecasts about the longterm evolution of the energy demand, resource commodities price and availability. It assumes that these forecasts are accurate and use this information to build a system deployment plan using a deterministic optimization approach. However, these long-term forecasts are subject to many uncertainties. For instance, the forecasted energy demand which, in our case studies, features a clear upward trend over the years might in fact decrease drastically, stay stable or simply increase in an smaller extent. In these cases, the LMES designed using the initial forecasts might be significantly oversized, leading to difficulties for meeting the demand in some time-steps (due to the minimum working load of the conversion devices) and to investment costs higher than necessary. This problem may be partly mitigated through the use of a stochastic optimization approach such as multi-stage stochastic programming. Multi-stage stochastic programming is a stochastic optimization paradigm based on the assumptions that a probabilistic description of the uncertain parameters is available and that some decisions do not have to be made right now but rather might be differed to a later point in time. When applied to our LMES design problem, this would mean that we have on hand some probabilistic information about the long-term demand forecast errors and that the design decisions relative to a future investment phase do not have to be made now but rather can be postponed to the beginning of this phase, at a point in time where more accurate information about the energy demands and resource commodities prices will be available. Such a multi-stage stochastic programming approach would typically rely on a scenario tree to represent the potential evolution of the input parameters over the investment phases. It would result in the formulation of an MILP displaying the same overall structure as the one investigated in this work. However, the size of this MILP will be drastically increased as it will be broadly proportional to the number of nodes in the scenario tree rather than to the number of investment phases involved in the investment planning horizon. Solving this MILP efficiently will thus be particularly challenging.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Research perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . Example of a local multi-energy system [Xie+20] . . . . . . . . . . . . . 3.1 RES diagram (above) and superstructure (bottom) of a simple cogeneration system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 System layout for a given year and operation schedule for a given hour for a simple cogeneration system . . . . . . . . . . . . . . . . . . . 3.3 RES of the studied district cooling systems . . . . . . . . . . . . . . . . 3.4 RES of the trigeneration system . . . . . . . . . . . . . . . . . . . . . . . 4.1 Representative days (2 typical days and 4 extreme days) for the cooling demand during Investment Phase 1 for the DCS located in City B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Piece-wise linear approximation of the performance curve at 30°C of a single-mode electric chiller . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Aggregate performance curves of single-mode electric chillers . . . . . 4.4 Inflow and outflow variables related to intermediate commodity c 1 = (n 1 , s 1 ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Example of inputs and outputs variables . . . . . . . . . . . . . . . . . 4.6 City A project: daily variations of the price for commodity (ELEC,1) . . 4.7 City A project: daily variations of the cooling demand during the 4 extreme days of Phase 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 City A project: performance curves of the single-mode electric chillers (also termed standard chillers) at 30 • C . . . . . . . . . . . . . . . . . . . 4.9 City A project: performance curves of the dual-mode electric chillers (also termed ice chillers) at 30 • C . . . . . . . . . . . . . . . . . . . . . . 4.10 City B project: daily variations of the price for commodity (ELEC, 1) . 4.11 City B project: daily variations of the cooling demand during the 4 extreme days of Phase 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 Trigeneration system: daily variations of the price for commodity (ELEC, 1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 Trigeneration: daily variations of the cooling and heating demand during the 4 corresponding extreme days of Phase 4 . . . . . . . . . . 4.14 Trigeneration system: performance curves of the single-mode electric chillers (also termed standard chillers) at 30 • C . . . . . . . . . . . . . . 4.15 Trigeneration system: performance curves of the dual-mode electric chillers (also termed ice chillers) at 30 • C . . . . . . . . . . . . . . . . . . 5.1 Example: Resolution of an integer program by a B&C algorithm with systematic rejection of the incumbent . . . . . . . . . . . . . . . . . . . . 5.2 Flow chart of hierarchical decomposition algorithm . . . . . . . . . . . 6.1 Operation cost of day (4,143) in the DCS of City B as a function of the number of installed chillers . . . . . . . . . . . . . . . . . . . . . . . . .

  constant Energy Type: S = solar radiation, E = electricity, G = natural gas, W = wind power, H = heating power, C = cooling power, ES = electricity sold to the grid, HW = hot water. Conversion Technology: CHP = combined heating and power, ST = solar thermal collector, EH = electric heater, CCHP = Combined cooling, heating and power, FC = fuel cell, HP = heat pump, WT = wind turbine, EC = electric chiller, EL = electrolyzer, Airco = air conditioning, GFAC = gas-fired absorption chiller. Storage Technology: H = heat storage, C = cooling storage, B = battery, H 2 = hydrogen.

FIGURE 3

 3 FIGURE 3.3: RES of the studied district cooling systems

FIGURE 3 . 4 :

 34 FIGURE 3.4: RES of the trigeneration system

FIGURE 4 . 1 :

 41 FIGURE 4.1: Representative days (2 typical days and 4 extreme days) for the cooling demand during Investment Phase 1 for the DCS located in City B

FIGURE 4 .

 4 FIGURE 4.2: Piece-wise linear approximation of the performance curve at 30°C of a single-mode electric chiller

  σ = P S . Let b be the index of the breakpoint of function π such that P σ ∈ [a b , a b+1 ]. The energy consumed by each device σ is thus given by: Q σ = s b P σ + c b where s b and c b are the slope and constant value of the b th line segment of π. The total energy consumed by the S devices is thus equal to Q = s b P + c b S. This equality holds for any value of P such that P S ∈ [a b , a b+1 ], i.e. any value of P ∈ [Sa b , Sa b+1 ]. This means that Π S is linear over the segment [Sa b , Sa b+1 ], with a slope equal to s b and a constant value of c b S.

FIGURE 4 .

 4 FIGURE 4.3: Aggregate performance curves of single-mode electric chillers

  step c . The largest number of discrete units of power that may be contracted for commodity c is thus MaxUnits c = ⌊ MaxC c C step c ⌋. We denote by SC ϕ,c the cost of contracting one discrete unit of size C step c of maximum nstantaneous power.

Figure 4 .

 4 4 illustrates this definition on a small RES diagram in which c 1 = (n 1 , s 1 ) is an intermediate commodity. The arrows in the diagram indicates that technologies p 1 and p 2 produce commodity c 1 and technology p 3 consumes commodity c 1 .

  Set of technologies consuming commodity c. P in f low cSet of technologies producing commodity c.Data EP t,cBuying (resp. selling) price for commodityc ∈ C R (resp. c ∈ C S,Sell ) in t ϵ t,HEAT Ratio of heat dissipation in t Dem t,c Demand for commodity c ∈ C S in t Conso max t,c Available amount of commodity c ∈ C R,Ph in t C step cSize of the discrete units used to define the value of the maximum instantaneous power contracted with the supplier of commodity c ∈ C R,Co SC ϕ,cCost of contracting one discrete unit of sizeC step c of commodity c ∈ C R,Co MaxUnits cMaximum number of discrete units of power that may contracted with the supplier of commodity c ∈ C R,Co Variables Cdis ϕ,c Number of discrete units of power contracted at phase ϕ with the supplier of commodity c ∈ C R,Co (Integer) Outflow t,c,p Amount of commodity c ∈ C consumed by technology p within t (Continuous) Inflow t,c,p Amount of commodity c ∈ C produced by technology p within t (Continuous) Conso t,c

FIGURE 4 . 4 :

 44 FIGURE 4.4: Inflow and outflow variables related to intermediate commodity c 1 = (n 1 , s 1 )

  c Inflow t,c,p = Dem t,c + Sold t,c ∀c ∈ C S , ∀t

  0 ∀c ∈ C S \ C S,Sell , ∀t (4.7) Constraints for the intermediate commodities An intermediate commodity c ∈ C I is produced and consumed within the system. The total amount produced by the technologies in P in f low cof production should thus be equal to the total amount consumed by the technologies in P

o

  t,COLD,m,b+1 -o t,COLD,m,b a t,,COLD,mb+1 -a t,COLD,m,b is the slope and cst t,COLD,m,b = o t,COLD,m,bs t,COLD,m,b a t,m,COLD,b the constant value of the line segment between breakpoints b and b + 1 in the performance curve of as a single active device of type m.

Figure 4 . 5

 45 Figure 4.5 illustrates these constraints for a technology denoted by p 1 . p 1 consumes energy of forms N in p 1 = {n 1 , n 2 } and produces energy in the form N out p = {n 3 }. Two types of devices denoted by l 1 and l 2 are available for this technology.

  p

FIGURE 4 . 5 :

 45 FIGURE 4.5: Example of inputs and outputs variables

  Figure 4.7, we show the cooling demand profile of the four extreme days in Phase 3. Each day is divided into H = 24 time steps, each one lasting ∆t = 1 hour. The maximum (resp. minimum non-zero) daily and hourly demand are 2.09 GWh and 93.8 MWh (resp. 423 kWh and 1.33 kWh).

FIGURE 4 .

 4 FIGURE 4.7: City A project: daily variations of the cooling demand during the 4 extreme days of Phase 3

  and 4.9. In order to build the piece-wise linear approximation of the performance curve to be used for each chiller of type m, each form of energy n and each timestep t, we use B t,n,m = 4 breakpoints. The breakpoint coordinates are obtained as described in Subsection 4.2.3.

FIGURE 4 .

 4 FIGURE 4.8: City A project: performance curves of the single-mode electric chillers (also termed standard chillers) at 30 • C

  Figure 3.3. Similar to the first case study, we have the set of available technologies P = {SMEC, DMEC, STO_COLD} and the set of energy forms N = {ELEC, COLD, ICE}. The set of resource commodities is C R = {(ELEC, 1)}, the set of intermediate commodities C I = {(ICE, 1)} and the set of supply commodities C S = {(COLD, 1)}.

FIGURE 4 .

 4 FIGURE 4.9: City A project: performance curves of the dual-mode electric chillers (also termed ice chillers) at 30 • C

FIGURE 4 .

 4 FIGURE 4.10: City B project: daily variations of the price for commodity (ELEC, 1)

FIGURE 4 .

 4 FIGURE 4.11: City B project: daily variations of the cooling demand during the 4 extreme days of Phase 5

FIGURE 4 .

 4 FIGURE 4.12: Trigeneration system: daily variations of the price for commodity (ELEC, 1)

FIGURE 4 .

 4 FIGURE 4.14: Trigeneration system: performance curves of the singlemode electric chillers (also termed standard chillers) at 30 • C

FIGURE 4 .

 4 FIGURE 4.15: Trigeneration system: performance curves of the dualmode electric chillers (also termed ice chillers) at 30 • C

Finally, the parameters

  relative to the energy storage technologies P STO = {STO_COLD, STO_HEAT} are the following.The capacity of the heat storage tank m = (STO_HEAT, 1) may be chosen at any value between 0kWh and MaxSize m = 180MWh. We use a discrete capacity unit Size step m equal to 9000kWh so that MaxNum m = 20. The cost for installing a kWh of heat storage capacity is 17CNY, which gives FC m = 17 × Size step m = 153kCNY. The capacity of the ice storage tank m = (STO_COLD, 1) may be chosen at any value between 0kWh and MaxSize m = 180MWh. We use a discrete capacity unit Size step m equal to 9000kWh so that MaxNum m = 20. The cost for installing a kWh of ice storage capacity is 120CNY, which gives FC m = 120 × Size step m = 1.08MCNY.

FIGURE 5 .

 5 FIGURE 5.1: Example: Resolution of an integer program by a B&C algorithm with systematic rejection of the incumbent

FIGURE 6 .

 6 FIGURE 6.1: Operation cost of day (4,143) in the DCS of City B as a function of the number of installed chillers

  as the optimal design decision and UB k as the minimum cost;

Algorithm 1 :

 1 Branch and Bound Algorithm for a minimization problem Data: The MILP formulation of the problem Result: The optimal solution and the minimum objective value 1 Solve the LP relaxation LP 0 of the original MILP problem; 2 if LP 0 is infeasible then

	3	return the original problem is infeasible;
	4 else
	5	if the solution of LP 0 is integer feasible then

6 return the

  solution of LP 0 as the optimal solution and the objective value of LP 0 as the minimum objective value ;Create the root node Node 0 of the search tree and add it to the list of active nodes;

	7	else
	8
	28	end
	29	else
	30	if Z n < Z best then

9 end end

Set Z best to positive infinity and x best to be empty; while the list of active nodes is not empty do Select an active node, Node i associated to problem LP i , remove Node i from the list of active nodes; Select an integer variable with a fractional value in the solution x i :

x j =

xj,i ; Create two child nodes of Node i in the search tree and index them by l i and r i ; Associate to Node l i an LP problem LP l i corresponding to problem LP i in which the constraint x j ≤ ⌊ xj,i ⌋ has been added to the formulation, and associate to Node r i an LP problem LP r i corresponding to problem LP i in which the constraint x j ≥ ⌈ xj,i ⌉ has been added to the formulation; for each newly created node n ∈ {l i , r i } do Solve the corresponding problem LP n ; if LP n is infeasible then 20 Prune Node n; else 22 Find the optimal solution x * n and the optimal objective value Z n ; 23 if x * n is not integer feasible then 24 if Z n ≥ Z best then 25 Prune Node n; 26 else 27 Add Node n to the list of active nodes; 31 Update the incumbent and set x best ← x * n and Z best ← Z n ; 32 end 33 end end end end return x best and Z best ;

  .2 partially describes (see the red numbers) this schedule by displaying the number of units of each type active during this hour together with the amount of thermal energy stored in the tank at the end of this hour.

	G A S	H E A T 1	H E A T 2	E L E C
	CHP
	Type1: 0/5
	Type2: 2/5
		STO
		Capacity: 30/50 MWh
	Boiler
	Type1: 8/10
	Type2: 2/10
	G A S	H E A T 1	H E A T 2	E L E C
	CHP
	Type1: 0/0/5
	Type2: 2/2/5
		STO
		Capacity: 25/30/50 MWh
	Boiler
	Type1: 6/8/10
	Type2: 0/2/10

FIGURE 3.2: System layout for a given year and operation schedule for a given hour for a simple cogeneration system

TABLE 4

 4

	.1: Notations for commodities

t,n,m,b

  Size of a discrete capacity unit for the device of type m = (p, l) belonging to technology p ∈ P CONT MaxNum m Maximum number of devices / discrete capacity units of type m that

	Set of indices (produced by the CCHP units in cooling mode) at time-step t is given by the follow-
	P ing relation.		Set of technologies
	N in p N out p C in n,p	Set of energy forms consumed by technology p ∈ P Set of energy forms produced by technology p ∈ P Set of commodities with energy form n consumed by technology p t,GAS,m ≥ P in P out t,HEAT,m η t,HEAT,m + P out t,COLD,m η t,COLD,m ∀m ∈ M CCHP , ∀t (4.14)
	C out n,p			Set of commodities with energy form n produced by technology p
	P DIS		Set of discrete technologies
	P CONT		Set of continuous technologies
	P CONV		Set of energy conversion technologies
	P STO		Set of energy storage technologies
	M			Set of device types
	M CONV		Set of energy conversion device types
	M STO		Set of energy storage device types
	M p			Set of device types belonging to technology p
	Data
	FC m		Total cost for installing one device / one discrete unit capacity of type
				m
	MC ϕ,m		Accumulated maintenance cost of one device / one discrete unit capac-
				ity of type m from phase ϕ to the end of system lifetime
	MaxSize m	Maximum size for the device of type m = (p, l) belonging to technol-
				ogy p ∈ P CONT
	Size	step m
				may be installed
	P min n,m		Minimum amount of energy of type n that a device of type m can pro-
				duce within one time-step
	P max n,m		Maximum amount of energy of type n that a device of type m can pro-
				duce within one time-step
	η t,n,m		Load-independent conversion efficiency of a device of type m ∈
				M Ordinate of breakpoint b of this piece-wise linear approximation
	Variables
	SD ϕ,m		Number of devices / number of discrete units of type m ∈ M newly
				added in the system at the beginning of phase ϕ (Integer)
	P in t,n,m		Total amount of energy of type n ∈ N in p consumed by the devices of
				type m ∈ M in t (Continuous)
	P out t,n,m		Total amount of energy of type n ∈ N out p	produced by the devices of
				type m ∈ M in t (Continuous)
	S t,m			Number of devices of type m ∈ M CONV active in t (Integer)
	Mode t,n,m	Number of active devices of type m ∈ M CONV operating under mode
				n ∈ N out p in t (Integer)
	STO t,m		Amount of energy stored in the devices of type m ∈ M STO at the be-
				ginning of t (Continuous)
				TABLE 4.2: Notations for technologies

t,GAS,m = P out t,ELEC,m η t,ELEC,m ∀m ∈ M CCHP , ∀t (4.13) Moreover, the relation between the amount of consumed gas and the amounts of heating energy (produced by the CCHP units in heating mode) and cooling energy CONV producing commodity n at time-step t B t,n,m Number of breakpoints in the piece-wise linear approximation of the performance curve for a device of type m ∈ M CONV producing commodity n with a load-dependent conversion efficiency at time-step t a t,n,m,b Abscissa of breakpoint b of this piece-wise linear approximation o

  Finally, for each representative day (ϕ, d) ∈ D, we have:

	STO ϕ,d,0,m = STO ϕ,d,H,m	∀m ∈ M STO , ∀(ϕ, d) ∈ D	(4.32)
	d,h,m + P in ϕ,d,h,n in p ,m -P out ϕ,d,h,n out p ,m = STO ϕ,d,h+1,m ∀m ∈ M STO , ∀(ϕ, d) ∈ D,
		∀h = 0, . . . , H -1	(4.31)

Constraints (4.31) state that the amount of energy stored at the beginning of hour h + 1, STO ϕ,d,h+1,m , is equal to the amount of energy already stored at the beginning of time-step h, STO ϕ,d,h,m , plus the total amount of energy put into the device during time-step h minus the energy released during h.

TABLE 4 .

 4 .3.

	Phase	1	2	3
	Length (years)	1	1	28
	Total Yearly Demand (GWh)	54.7 250.0 276.8
	Max Hourly Demand (MWh) 13.6 62.2	93.8

3: City A project: upward trend of the annual demand

TABLE 4 .

 4 4: City A project: available types of chiller SMEC and DMEC at an ambient temperature of 30 • C are displayed in Figures 4.8

	m (SMEC,1)	P max COLD,m (MW) P max ICE,m (MW) FC m (CNY) 9 -1.3 × 10 8
	(SMEC,2)	9	-	1.2 × 10 8
	(SMEC,3)	5	-	7.9 × 10 7
	(DMEC,1)	8	6	1.3 × 10 8
	(DMEC,2)	5	3	8.2 × 10 7

max n,m . The performance curves of the available types of

Table 4 .

 4 5.

	Phase	1	2	3	4	5
	Length(yr)	1	1	1	1	26
	Total Yearly Demand(GWh)	0.8 8.4	26.1 48.8 62.2
	Max Hourly Demand(MWh) 1.4 14.0 45.3 90.0 132.0

TABLE 4 .

 4 5: City B project: upward trend of the annual demand

TABLE 4 .

 4 • C are thus displayed in Figures 4.8 and 4.9. However, the installation cost FC m in City B differs from the 6: City B project: available types of chiller one in City A project and is provided in Table 4.6.

	(SMEC,1)	P max COLD,m (MW) P max ICE,m (MW) FC m (CNY) 9 -1.3 × 10 8
	(SMEC,2)	9	-	1.2 × 10 8
	(DMEC,1)	8	6	1.3 × 10 8

TABLE 4 .

 4

	Phase	1	2	3	4
	Length (years)	1	1	1	27
	Total Yearly Cooling Demand (GWh)	2.3	3.1	12.7 18.9
	Max Hourly Cooling Demand (MWh) 4.6	5.8	20.2 30.2
	Total Yearly Heating Demand (GWh)	14.0 24.2 31.1 38.3
	Max Hourly Heating Demand (MWh) 11.0 17.3 24.6 32.6

7: Trigeneration system: upward trend of the annual demand

  * PC m ∑ Φ φ=1 β φ .

	m (CCHP,1)	P max COLD,m 1.5	P max ICE,m -	P max HEAT,m 1.5	P max ELEC,m 1.5	FC m (CNY) 1.70 × 10 7
	(BOILER,1)	-	-	7	-	2.94 × 10 6
	(SMEC,1)	3.5	-	-	-	5.20 × 10 6
	(SMEC,2)	3.5	-	-	-	3.05 × 10 6
	(DMEC,1)	3.5	2.4	-	-	4.39 × 10 6
	(ASHP,1)	1.6	-	1.6	-	4.38 × 10 6

TABLE 4

 4

.8: Trigeneration system: available types of technology

TABLE 4 .

 4

		10: City B project: direct resolution by CPLEX 12.8 solver
		|D ϕ |	12	16	20	24	28	32	36
	MILP	#Var	76080 101424 123600 148944 175872 201216 226560
	size	#IntVar	16164	21540	26244	31620	37332	42708	48084
		#Cons	86420 115220 140420 169220 199820 228620 257420
	Direct	Gap LP	5.11% 5.17%	5.20%	5.21%	5.15%	5.14%	5.13%
	resolution Time	5934s * 6302s *	7200s 5259s * 6155s * 5968s * 5327s *
		Gap MIP 0.06% 0.26%	0.23%	0.35%	0.32%	0.16%	1.21%

"*" means that the computer ran out of memory before reaching the time limit.

TABLE 4

 4

.11: Trigeneration project: direct resolution by CPLEX 12.8 solver

  ). Let Z best CP denote the value of the best feasible solution of CP (x best , y best , z best ) found since the beginning of the algorithm, which is also called the incumbent of CP. If Z ♯ CP and add Z best CP as a new upper cutoff value for SRP in the upper-level B&C algorithm SRP, i.e. add the constraint f 0

CP is lower than Z best CP , we set (x best , y best , z best ) to (x ♯ , y ♯ , z ♯ ) and Z best CP to Z ♯

  .., x 2 Φ ) is a potential design found after x 1 during the upper level branch-and-cut search. For each phase ϕ and each d ∈ D ϕ , we check whether (x 2 ϕ , ϕ, d) is a key of the hash table. If so, i.e. if x 2 ϕ = x 1 ϕ , the corresponding relaxed operation cost f ϕ,d ( ỹϕ,d , zϕ,d ) can be directly replaced in the computation of Z 2 CP by the value Z ϕ,d (x 1 ϕ ) stored in the hash table so that we avoid a new resolution of OP ϕ,d (x 1 ϕ ) = OP ϕ,d (x 2 ϕ ). If (x 2 ϕ , ϕ, d) is not a key of the hash table, we solve operation sub-problem OP ϕ,d (x 2 ϕ

  ♯ , y ♯ , z ♯ ) to (x best , y best , z best ) and Z ♯ best , y best , z best ) ← (x ♯ , y ♯ , z ♯ ), set Z best ← Z ♯

		Start	Set D tmp ← D
	S,Sell , add the valid inequalities (5.29)-(5.30) to the formulation of SRP; 2 Set Z best CP to positive infinity; 3 Create the root node of the B&C search tree ; 4 while active nodes exist in B&C search tree do 5 Branch from an active node in the tree ; 6 if an integer feasible solution (x ♯ , ỹ♯ , z♯ ) with an objective value Z ♯ SRP ≤ Z best CP is found then 7 Set IsNewIncumbent to True; 8 for each phase ϕ ∈ {1, ..., Φ} and each day d ∈ D ϕ do 9 if (x ♯ ϕ , ϕ, d) is in the hash table then 10 Retrieve the corresponding value of Z ϕ,d (x ♯ ϕ ); 11 Set Z ♯ CP to Z ♯ CP + Z ϕ,d (x ♯ ϕ ) -f ϕ,d ( ỹ♯ ϕ,d , z♯ ϕ,d ); 12 else 13 Solve Problem OP ϕ,d (x ♯ ϕ ) as an MILP; 14 if OP ϕ,d (x ♯ ϕ ) is infeasible then 15 Add the no-good cut (5.28) corresponding to x ♯ ϕ to SRP; 17 break ; /* To break the for loop in Line 8 */ 18 else Add valid inequal-ities (5.29)-(5.30) to formulation SRP and set Z best ← +∞ Carry out a B&B search based on formulation SRP using Z best as cutoff value Any active node in the B&B tree? End Select an active node as current node Current solution inte-ger feasible for SRP? Obtain (x ♯ , ỹ♯ , z♯ ) and set Z ♯ CP ← Z ♯ SRP D tmp = ∅? Select (ϕ, d) ∈ D tmp and set D tmp ← D tmp \{(ϕ, d)} (x ♯ , ϕ, d) in hash table? Retrieve Z ϕ,d (x ♯ ϕ ) from hash table Z ♯ CP ← Z ♯ CP + Z ϕ,d (x ♯ ϕ ) -f ϕ,d ( ỹ♯ , z♯ ) Z ♯ CP < Z best ? Reject current solution and add the no-good cut coressponding to x ♯ to formulation SRP Solve the subprblem OP ϕ,d (x ♯ ϕ ) with for-mulation (5.16)-(5.20) Is OP ϕ,d (x ♯ ϕ ) feasible? Store Z ϕ,d (x ♯ ) in hash table with key (x ♯ ϕ , ϕ, d) Add the no-good cut corresponding to x ♯ ϕ to formulation SRP Set (x FIGURE 5.2: Flow chart of hierarchical decomposition algorithm
	19	Store the value of Z ϕ,d (x ♯ ϕ ) under the new key (x ♯ ϕ , ϕ, d) in
		the hash table;
	20	Set Z ♯ CP to Z ♯ CP + Z ϕ,d (x ♯ ϕ ) -f ϕ,d ( ỹ♯ ϕ,d , z♯ ϕ,d );
	21	end
	22	end
	23	if Z ♯ CP > Z best CP then
	24	Set IsNewIncumbent to False;
	26	Add the no-good cut (5.27);
	27	break;	/* To break the for loop in Line 8 */
	28	end
	29	end
	30	if IsNewIncumbent is True then
	31 32	Set (x CP to Z best CP ; Add cut f 0 (x) + ∑ (ϕ,d)∈D f ϕ,d (y ϕ,d , z ϕ,d ) ≤ Z best CP to SRP ;
	33	Reject (x ♯ , ỹ♯ , z♯ );
	34	Add the no-good cut (5.27);
	35	end
	36	end
	37 end
	38 return (x best , y best , z best ) and Z best CP ;

16

Set IsNewIncumbent to False; 25 Reject (x ♯ , ỹ♯ , z♯ );

TABLE 5

 5

		.1: City A project: numerical comparison of the solution
			methods
		|D ϕ |	6	8	10	12	14
	MILP	#Var	13820 18380 22940 27500	32060
	size	#IntVar	3625	4825	6025	7225	8425
		#Cons	23084 30764 38444 46124	53804
	Direct	Gap LP	5.32% 4.26% 4.58% 4.77% 4.38%
	resolution Time	7200s 7200s 7200s 2589s * 4411s *
		Gap MIP	0.01% 0.02% 0.02% 0.04% 0.06%
	Initial	Gap LP	5.32% 4.26% 4.58% 4.77% 4.38%
	Hier. Dec. #IFS SRP	4087	1603	1303	2348	1688
		#OP solved 45406 35193 31752 55160	34656
		Time	5162s 7200s 7200s 7200s	7200s
		Gap MIP	0.00% 0.92% 0.86% 1.60% 0.66%
	Extended Gap LP	4.13% 3.37% 4.10% 4.24% 4.04%
	Hier. Dec. #IFS SRP	1764	1432	1441	394	1660
		#OP solved	505	794	618	520	1287
		Time	954s	2889s 4580s 7200s	4062s
		Gap MIP	0.00% 0.00% 0.00% 0.07% 0.00%
	"*" means that the computer ran out of memory before reaching the time limit.

TABLE 5

 5

	.2: City B project: numerical comparison of the solution
	methods

  . Two system layouts are said to be comparable if one is smaller than the other.Moreover, we note that the constraints (6.3) coupling the design and the operation variables in CP have a specific structure. Namely, they either take the form y ϕ,d,h,j ≤ x ϕ,i (see e.g. Constraints (4.9),(4.15), (4.19), (4.22), (4.25) for the energy conversion technologies) or the form z ϕ,d,h,j ≤ c i x ϕ,i with c i > 0 (see e.g. Constraints (4.5) for the resource commodities in C R,Co and Constraints (4.29) for the energy storage technologies.) Let us now consider two operation sub-problems OP ϕ,d (x ♯ ϕ ) and OP ϕ,d (x • ϕ ) and their feasible spaces denoted by V ♯ ϕ,d and V • ϕ,d . Thanks to the key observation on the coupling constraints mentioned above, we note that if x ♯ ϕ,i ⪯ x • ϕ,i , all the coupling constraints involved in OP ϕ,d (x ♯ ϕ ) have a right-hand side smaller than the one of the coupling constraints involved in OP ϕ,d (x • ϕ ) and are thus stricter so that V ♯ ϕ,d ⊆ V • ϕ,d . This inclusion relationship between their feasible regions implies that if both operation scheduling sub-problems OP ϕ,d (x ♯ ϕ ) and OP ϕ,d (x • ϕ ) are feasible with optimal objective value Z ♯

• ϕ if x ♯ ϕ is element-wise smaller than x • ϕ , i.e. if we have x ♯ ϕ,i ≤ x • ϕ,i for each i = 1 . . . ν ϕ . In this case, we note x ♯ ϕ ⪯ x • ϕ ϕ,d and Z

  4,143 (x ♯ 4 ) is feasible for a given design vector x ♯ 4 , Z 4,143 (x ♯ 4 ) is a lower bound for the optimal value of OP 4,143 (x 4 ) as long as x 4 ⪯ x ♯ 4 . See e.g. how Z 4,143 (10, 10) = 357996 is a lower bound for Z 4,143 (x 4 ) as long as x 4 ⪯ (10, 10).

	Let x ♯

ϕ ∈ proj x ϕ (F ϕ ) be a system layout at phase ϕ feasible for each day d ∈ D ϕ .

19 end 20 end 21 if IsFeasible is False then 22 break

	15	else
	16	Add the feasibility cut (6.28) corresponding to x k ϕ to MP k ;
		;	/* To break the for loop started in Line 9 */
	23	else
	24	Add the optimality cut (6.27) corresponding to x k ϕ to MP k ;
	25	end
	26	end
	27	if IsFeasible is True then
	28

d (x k ϕ ); 17 Set IsFeasible to False; 18 break ; /* To break the for loop started in Line 11 */

TABLE 6 .

 6

1: City A project: numerical comparison of the solution methods

TABLE 6

 6

		.2: City B project: numerical comparison of the solution
				methods
		|D ϕ |	22	30	38	58	78	165
	Generalized #Iter	4	4	4	4	4	4
	Benders	Time	207s	364s	463s	991s	1593s 4370s
		Gap MIP 0.00% 0.00% 0.00% 0.00% 0.00% 0.00%
	"*" includes all days with non-zero cooling demand in one year

TABLE 6 .

 6 4127s * 2929s * 3244s * 3863s * 5060s * 3713s *

	3: City B project: results obtained with the generalized Ben-
	ders' decomposition algorithm on large instances

*

"*" means that the computer ran out of memory before reaching the time limit.

TABLE 6 .

 6

4: Trigeneration project: numerical comparison of the solution methods

TABLE 6 .

 6 7: City B project: deployment plan obtained as a function of the number of representative days obtain, for each phase ϕ, a total 'actual' operation cost Z actual m∈M ∑ c∈C R EP ϕ,d,h,c Conso actual ϕ,d,h,c . The corresponding results are displayed in Table

	|D ϕ |	6	14	22	30	38	58	78
	Obj. value	5.50 × 10 8	5.74 × 10 8	5.74 × 10 8	5.75 × 10 8	5.74 × 10 8	5.73 × 10 8	5.74 × 10 8
	Design cost	2.44 × 10 8	2.36 × 10 8	2.36 × 10 8	2.36 × 10 8	2.36 × 10 8	2.36 × 10 8	2.36 × 10 8
	Operation cost	3.06 × 10 8	3.38 × 10 8	3.39 × 10 8	3.40 × 10 8	3.39 × 10 8	3.37 × 10 8	3.38 × 10 8

TABLE 6

 6 × 10 6 -0.75% 4.38 × 10 6 -0.83% 4.41 × 6 2 1.89 × 10 7 -0.90% 1.88 × 10 7 -1.12% 1.90 × 7 3 3.15 × 10 8 -0.70% 3.14 × 10 8 -0.87% 3.17 × 8 Total 3.38 × 10 8 -0.71% 3.38 × 10 8 -0.88% 3.41 × 8 × 10 4 -0.41% 4.98 × 4 2 5.85 × 10 5 -0.86% 5.88 × 10 5 -0.26% 5.90 × 5 3 1.63 × 10 6 -0.22% 1.62 × 10 6 -0.38% 1.63 × 6 4 2.82 × 10 6 -0.75% 2.83 × 10 6 -0.45% 2.84 × 6 5 3.76 × 10 7 -0.98% 3.77 × 10 7 -0.68% 3.79 × 7 Total 4.26 × 10 7 -0.93% 4.28 × 10 7 -0.65% 4.30 × 7

	Phase Z ϕ optim	ϵ optim	Z sim ϕ	ϵ sim	Z actual ϕ
	1	4.38 .9: City A project: comparison of the approximated, simu-
			lated and actual operation cost
	Phase Z ϕ optim	ϵ optim	Z sim ϕ	ϵ sim	Z actual ϕ
	1	5.17 × 10 4 3.85%	4.96

TABLE 6 .

 6 10: City B project: comparison of the approximated, simulated and actual operation cost

FIGURE 4.6: City A project: daily variations of the price for commodity (ELEC,1)

per contracted kW for is equal to 270CNY for all phase ϕ, which gives SC (ELEC,1),ϕ = 270 × C Let us now consider the available conversion technologies P CONV = {CCHP, BOILER, ASHP, SMEC, DMEC}.

A single type of CCHP unit (L CCHP = 1) is available for installation in the LMES. The available model has an electricity generation capacity equal to P max ELEC,(CCHP,1) = 1500kW and thermal generation capacities equal to P max COLD,(CCHP,1) = P max HEAT,(CCHP,1) = 1500kW. Its minimum load rate is 50%, which means that an active CCHP unit should produce at least P min ELEC,(CCHP,1) = P min COLD,(CCHP,1) = P min HEAT,(CCHP,1) = 750kW of electric power. Its electric and thermal conversion efficiencies are load-independent and constant: η t,ELEC,(CCHP,1) = η t,COLD,(CCHP,1) = η t,HEAT,(CCHP,1) = 38% for all t.

The purchasing price of one unit of (CCHP, 1) is 7MCNY and the corresponding total installation cost is FC (CCHP,1) = 17MCNY.

There is also a single type of boiler: L CCHP = 1. The available model may be operated at any value between P min HEAT,(BOILER,1) = 0kW and its maximum capacity P max HEAT,(BOILER,1) = 7000kW. Its conversion efficiency coefficient is load-independent and constant: η t,HEAT,(BOILER,1) = 95% for all t. The cost for purchasing a unit of (BOILER, 1) is 1.1MCNY and the total installation cost per unit is FC (BOILER,1) = 2.94MCNY.

As for the single-mode electric chillers, L SMEC = 2 types of devices are available. Both types of chiller have a maximum working load of P max COLD,(SMEC,1) = P max COLD,(SMEC,2) = 3.5MW. Chillers of type (SMEC, 1) have a minimum load rate of 10% so that P min COLD,(SMEC,1) = 0.35MW. Chillers of type (SMEC, 2) have a minimum load rate of 20% so that P min COLD,(SMEC,2) = 0.70MW. Their conversion efficiencies are load-and temperature dependant. They are depicted by a set of convex performance curves providing the consumed electric power as a function of the produced cooling power under different ambient temperatures. The performance curves of both types of chillers at an ambient temperature of 30 • C are shown in Figure 4.14. The purchase price of one device of type (SMEC, 1) is 2.1MCNY and the total installation cost is FC (SMEC,1) = 5.2MCNY. The purchase price of one device of type (SMEC, 1) is 1.2MCNY and the total installation cost is FC (SMEC,2) = 3.1MCNY.

A single type of dual-mode chiller is available for installation in the system: L DMEC = 1. Its maximum cooling capacity is P max COLD,(DMEC,1) = 3.5MW and its maximum ice-making capacity is P max ICE,(DMEC,1) = 2.4MW. The minimum load rate for Chapter 6 Generalized Benders' decomposition algorithm

Introduction

The numerical results presented in Chapter 5 showed that the extended hierarchical decomposition algorithm is more effective at solving the LMES optimal design problem than an off-the-shelf MILP solver for the instances related to the DCS projects. However, numerical difficulties were still encountered while trying to solve the instances related to the trigeneration system. We thus investigate in this chapter another solution approach which is also based on a decomposition strategy. This approach can be seen as a generalized Benders' decomposition algorithm, more precisely as an extension of the conventional Benders' decomposition algorithm to a case where the second-stage sub-problems involve integer variables. In particular, we propose a new type of optimality and feasibility cuts exploiting the special structure of the constraints linking the design and operation variables in our problem.

Thanks to these new cuts, we are able to develop a generalized Benders' decomposition algorithm converging in a finite number of iterations to an optimal solution of the original problem CP.

In this chapter, we first present in Section 6.2 the overall framework used to develop a generalized Benders' decomposition algorithm for our problem. The optimal and finite convergence of this algorithm relies on the use of strong dual bounding functions. We explain in Section 6.3 how such functions can be obtained for our problem by exploiting the special structure of the constraints coupling the design and operation variables. The resulting master problem contains a set of non-linear constraints. Section 6.4 thus investigates how this master problem may be reformulated as an MILP thanks to the introduction of a set of additional binary variables. Section 6.5 summarizes the outline of the proposed algorithm. Finally, Section 6.4 reports the results of the computational experiments carried out to assess the numerical performance of the algorithm. It also presents the outcome of a postoptimization simulation study carried out to evaluate the quality of the deployment plans obtained for the two DCS projects.

A generalized Benders' decomposition framework

Bolusani and Ralphs [START_REF] Bolusani | A Framework for Generalized Benders' Decomposition and Its Application to Multilevel Optimization[END_REF] recently described a framework to extend the wellknown decomposition framework introduced by Benders [START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF] to a more general class of optimization problems. In particular, they highlighted the fact that Benders' decomposition may be used even if the second-stage sub-problems do not display Algorithm 3: Basic generalized Benders' decomposition algorithm Data: Parameters of commodities and technologies, construction phases, discount rate and a value of the optimality tolerance gap ϵ Result: The optimal system deployment plan and the minimum total cost 1 Set the number of iteration k to 0; 2 Set the upper bound UB 0 to +∞ and the lower bound LB 0 to -∞;

Solve the master problem MP k using formulation (6.25)-(6.34) with an MILP solver; than 12 days. Overall, these results indicate that, provided the number of representative days used to estimate the operation cost is large enough and the generalized Benders' decomposition algorithm can converge within the time limit, the design decisions do not depend on the cardinality of the subset of representative days and may thus be recommended for implementation.

We then seek to further assess the impact of using a limited number of representative days and a piece-wise linear approximation of the devices' performance curves to estimate the future operation cost in our optimization model. This assessment relies on a post-optimization numerical simulation study aiming at computing as accurately as possible the actual operation cost corresponding to the deployment plan recommended by our optimization model. Due to time reasons, we carried out this assessment only on the DCS projects. Note that in these projects, all the available conversion devices have non-linear convex performance curves.

For this simulation study, we use the following procedure:

1. We retrieve the design decisions for each phase ϕ obtained when solving the MILP problem formulated in Chapter 4 to optimality. For City A project (resp.

9 9 0 18 0 0 18 18 0 0 18 18 0 0 18 18 0 0 18 18 0 0 18 18 0 0 18 18 Heat Stor. cap. (MWh) 18 0 0 45 18 0 0 36 18 0 0 36 18 0 0 36 18 0 0 36 18 0 0 36 18 0 0 36 Cont. power (MW) 1 1 1.5 2 1 1 1.5 2 1 1 1.5 2 1 1 1.5 2 1 1 1.5 2 1 1 1.5 2 1 1 1.5 2 computed by the model for day d. We also record the leaving ice inventory level, i.e. STO sim ϕ,d,H-1,m . • When the scheduling problem of Day 364 and Day 365 is solved, which is the last two-day scheduling problem of the phase, the load allocation of each type of chiller of both two days, P out,sim ϕ,364,h,n,m and P out,sim ϕ,365,h,n,m , as well as the corresponding total operation cost, Z sim ϕ,364 and Z sim ϕ,365 , will be recorded.

We finally compute, for each phase ϕ, a total 'simulated' operation cost Z sim ϕ equal to Z sim ϕ = ∑ 365 d=1 Z sim ϕ,d .

3. We use the performance curves provided by the chillers' manufacturer to compute the actual operation cost. Thus, for each phase ϕ = 1, ..., Φ, each day d = 1, ..., 365, each hour h = 1, .., H, each type of chiller m ∈ M SMEC ∪ M DMEC and each operating mode n, we retrieve the value P out,sim ϕ,d,h,n,m obtained in Step 2 of the simulation procedure and use the chiller performance curve provided by its manufacturer to compute the corresponding actual electricity consumption P in,actual ϕ,d,h,ELEC,m . This allows us to compute the total consumption of commodity (ELEC, 1)for the corresponding time-step: Conso actual ϕ,d,h,(ELEC,1) . We thus