UNIVERSITÉ DE MONTRÉAL

OPTIMISATION STOCHASTIQUE DE LA RÉPARTITION SPATIO-TEMPORELLE
D'UN VOLUME D'EAU AUX GROUPES TURBO-ALTERNATEURS D'UN SYSTÈME
DE PRODUCTION HYDROÉLECTRIQUE

SARA SÉGUIN

DÉPARTEMENT DE MATHÉMATIQUES ET DE GÉNIE INDUSTRIEL
ÉCOLE POLYTECHNIQUE DE MONTRÉAL

THÈSE PRÉSENTÉE EN VUE DE L'OBTENTION
DU DIPLÔME DE PHILOSOPHIÆ DOCTOR
(MATHÉMATIQUES DE L'INGÉNIEUR)

AOÛT 2016
UNIVERSITÉ DE MONTRÉAL

ÉCOLE POLYTECHNIQUE DE MONTRÉAL

Cette thèse intitulée :

OPTIMISATION STOCHASTIQUE DE LA RÉPARTITION SPATIO-TEMPORELLE
D'UN VOLUME D'EAU AUX GROUPES TURBO-ALTERNATEURS D'UN SYSTÈME
DE PRODUCTION HYDROÉLECTRIQUE

présentée par : **SÉGUIN Sara**

en vue de l'obtention du diplôme de : **Philosophiæ Doctor**

a été dûment acceptée par le jury d'examen constitué de :

M. GENDREAU Michel, Ph. D., président

M. AUDET Charles, Ph. D., membre et directeur de recherche

M. CÔTÉ Pascal, Ph. D., membre et codirecteur de recherche

M. DESAULNIERS Guy, Ph. D., membre

M. SHAWWASH Ziad, Ph. D., membre externe
# DÉDICACE

À maman
et Nana
REMERCIEMENTS

Quelle aventure ! Je n'ai pas vu passer les quatres dernières années de ma vie !

J'ai eu une opportunité incroyable de trouver ce projet de recherche alors que j'étais blasée
au travail... quelques courriels, une rencontre, et j'étais lancée pour une maîtrise. Un an plus
tard, j'étais lancée pour un doctorat, et maintenant je suis lancée pour un post-doctorat...

Tout d'abord je tiens à remercier les deux personnes les plus importantes dans ce projet, mes
directeurs. Charles et Pascal, je ne saurais trouver les mots justes pour dire à quel point j'ai
apprécié travailler avec vous. Nous avons formé une belle équipe ultra efficace et je veux vous
en remercier. Lorsqu'on a du support et des encouragements de ses directeurs, je crois que
cela fait toute la différence et nous aide à continuer. Je ne me suis jamais sentie seule et vous
avez toujours été disponible pour répondre à mes question concernant le projet, mais aussi
mes question existentielles revenues quelques fois au fil des années. J'ai énormément appris de
votre expérience et de vos compétences et vous m'avez aidée à devenir une doctorante outillée
pour poursuivre mon chemin en recherche. Vous m'avez toujours encouragée à saisir les
opportunités qui se présentent et cette liberté aussi compte pour beaucoup lors d'un doctorat.
Bref, j'ai adoré faire ce projet avec vous et je recommencerais sans aucune hésitation demain
matin ! Je souhaite à tous les étudiants de trouver des directeurs comme vous (sérieusement !).

Ce projet a été réalisé en collaboration avec Rio Tinto au Saguenay où j'ai dû me déplacer
régulièrement. Merci à l'équipe de gestion des ressources hydriques pour votre bel accueil
et votre implication dans le projet. De plus, ce projet m'a fait découvrir la belle région du
Saguenay Lac-St-Jean. J'ai maintenant un beau réseau d'amis dans cette région et tout ça
grâce à Catherine, à Hubert et Léo aussi. Je me sens chez moi quand je suis chez toi Cath et
c'est toujours le fun ! Et à Stéphane pour le vélo et la raquette, et les bières !

J'ai eu la chance de passer 6 mois à Trondheim, en Norvège, afin de travailler sur mon 2e
article avec Stein-Erik Fleten. Ce fut une expérience intéressante autant sur le plan acadé-
mique/professionnel que sur le plan personnel. Travailler dans un autre pays permet d'élargir
ses horizons et son réseau de contact et cela a été très bénéfique pour la poursuite de mon
doctorat. J'ai rencontré des gens supers, autant à l'université qu'à l'extérieur. Tusen takk
à mon colocataire norvégien Bjørnar, Eve, Marianne, Steinar et Norvald pour l'immersion
complète dans ce beau pays.

Entre tous ces déplacements, il y a aussi mes collègues du GERAD qui ont été plus qu'im-
portants dans cette aventure. À mes amis qui comprennent que je parte et que je revienne
sans trop savoir quand, vous êtes toujours là pour mes folies et je vous en remercie !! En plus vous venez me visiter, que ce soit à Québec, au Saguenay ou en Norvège :)

Si j'écris ces remerciements aujourd'hui c'est grâce à ma mère. Je ne sais pas par où commencer! Tu m'as toujours soutenue dans tous mes projets et tu m'as toujours encouragée à continuer. Tu ne connais pas l'adversité et je n'ai jamais de ma vie rencontré quelqu'un d' aussi fort que toi. Tu es toujours souriante et prête à aider, merci pour absolument tout, pour avoir fait de moi la personne que je suis aujourd'hui. Bien que je ne vous voie pas souvent, merci à toute ma famille pour votre éternel support et votre compréhension.

Un doctorat c'est un accomplissement énorme, autant sur le plan personnel que professionnel, alors merci à tous ceux qui ont fait de cette aventure une expérience enrichissante.

Alors sur une note finale : yaya¹⁰⁰⁰.
RÉSUMÉ

Le procédé de fabrication de l'aluminium consomme énormément d'énergie. Le Canada est le 4e producteur mondial d'aluminium et 90% de la production est effectuée au Québec, notamment pour son abondance en ressources hydriques. Cette province produit de l'hydroélectricité à des prix compétitifs et plusieurs producteurs d'aluminium possèdent des installations au Québec. L'entreprise minière Rio Tinto opère des alumineries dans la région du Saguenay Lac-St-Jean et est aussi propriétaire d'un système de production hydroélectrique, lui permettant de satisfaire 90% de la demande en énergie de ses alumineries. L'entreprise achète la portion restante à Hydro-Québec et il est dans leur intérêt de gérer efficacement leur réseau hydroélectrique afin de produire le maximum d'énergie avec l'eau disponible.

Plusieurs modèles d'optimisation sont utilisés pour gérer efficacement un système de production hydroélectrique. L'un d'entre eux est le modèle à court terme, aussi nommé le problème de chargement et de répartition optimal. Un système de production hydroélectrique est composé de centrales, qui abritent les turbines. Il existe deux types de centrales : au fil de l'eau et à réservoir. Les centrales au fil de l'eau possèdent peu de réserve d'eau et la puissance générée par les turbines dépend principalement du débit de l'eau. Les centrales à réservoir sont situées à côté d'un barrage permettant de moduler le volume d'eau dans les réservoirs et ainsi la production de puissance. La solution obtenue par cette optimisation détermine le volume des réservoirs, les turbines en marche dans chaque centrale et le débit turbiné par chacune d'entre elle, pour chaque pas de temps, habituellement horaire ou journalier. Ce problème est complexe à résoudre, car les fonctions de production sont non convexes et non linéaires. De plus, chaque turbine a une efficacité différente et les démarrages des turbines doivent être limités puisqu'ils causent une usure prématurée de l'équipement. Enfin, les précipitations qui tombent sur les bassins hydrographiques sont incertaines, ce qui signifie que les apports naturels d'eau dans les réservoirs sont inconnus au moment de prendre une décision. Cette thèse propose une modélisation du modèle à court terme pour le réseau hydroélectrique du Saguenay Lac-St-Jean, ainsi qu'une méthode d'optimisation stochastique pour le résoudre.

Tout d'abord, une modélisation déterministe du problème est proposée, c'est-à-dire que les apports naturels sont considérés connus. Un algorithme de programmation dynamique permet de calculer des surfaces de puissance pour chaque centrale qui sont ensuite utilisées comme fonction objectif dans la méthode d'optimisation en deux phases développée. La première phase, le problème de chargement, résout un modèle non linéaire mixte en nombres entiers. La deuxième phase, le problème de répartition, résout un modèle linéaire en nombres entiers.
Les résultats obtenus avec la méthode proposée sont comparés aux décisions opérationnelles historiques et démontrent que la méthode permet d'augmenter la production d'énergie. De plus, une comparaison avec une modélisation en une seule phase confirme que la méthodologie proposée permet de résoudre rapidement le modèle à court terme.

Ensuite, la modélisation déterministe est étendue afin de considérer des apports naturels in-
certains et un modèle à court terme stochastique est développé. Une méthode de génération
d'arbres de scénarios, basée sur la minimisation de la distance imbriquée est utilisée afin de
créer les arbres de scénarios utilisés dans le modèle d'optimisation, à partir de prévisions
d'apports. L'équivalent déterministe du modèle stochastique est résolu, c'est-à-dire qu'un
modèle déterministe est résolu à chaque noeud de l'arbre de scénarios. Un horizon roulant
est retenu pour valider le modèle d'optimisation. Les prévisions d'apports sont mises à jour
quotidiennement, un arbre de scénarios est généré, le modèle stochastique est résolu, la so-
lution du premier noeud est implémentée, les volumes des réservoirs sont mis à jour avec
la vraie réalisation des apports et le processus est relancé pour tout l'horizon de planifica-
tion. Les résultats numériques comparent la solution obtenue par l'optimisation stochastique
à l'utilisation du scénario médian des prévisions dans l'horizon roulant et démontrent que
l'utilisation d'un modèle stochastique permet de produire plus d'énergie. Des tests avec des
arbres de scénarios comportant un nombre différent d'étapes et de scénarios sont effectués.

Finalement, la complexité des arbres de scénarios, utilisés dans le modèle d'optimisation à
court terme, est étudiée. La méthode de génération d'arbres de scénarios requiert des para-
mètres qui sont : le nombre d'étapes, le nombre de noeuds fils par étape et l'agrégation de
chaque étape. L'horizon roulant, comportant la génération des arbres de scénarios et l'opti-
misation, est défini comme un problème d'optimisation boîte noire et les paramètres de la
méthode de génération d'arbres de scénarios sont optimisés. Le problème consiste à maxi-
miser la production d'énergie de l'horizon roulant en considérant la valeur de l'eau restant
dans les réservoirs à la fin de l'horizon de planification et en pénalisant les démarrages de
turbines. La solution obtenue avec l'arbre permettant de maximiser l'énergie est comparée
avec l'optimisation utilisant des peignes de scénarios pour représenter l'incertitude des ap-
ports naturels. Les résultats numériques démontrent que, dans ce cas précis, l'utilisation d'un
éventail de peignes de scénarios permet d'obtenir des résultats comparables à un arbre de
scénarios complexe, mais que le temps de calcul est grandement diminué.
ABSTRACT

Aluminium production requires huge amounts of energy. Canada is the fourth largest alu-
minium producer in the world and 90% of the production takes place in Québec, mainly for
its abundant water resources. This province produces hydropower at competitive prices and
many aluminium producers operate in Québec. The mining company Rio Tinto is one of
them and owns many aluminium plants in the Saguenay Lac-St-Jean region. This company
also owns an hydroelectric production system, which allows them to satisfy the demand in
energy of the aluminium plants at a 90% fulfillment. For the remainder, they need to buy
the energy from Hydro-Québec so it is at their advantage to operate the hydropower system
as efficiently as possible, in order to minimize their energy purchase.

Many optimization models are used to manage efficiently hydropower systems. One of them is the short-term, or unit commitment and loading problem. An hydropower system is composed of power plants which contains the turbines. There are two types of power plants: run of the river and power plants with storage reservoir. Run of the river plants have small water reserve capacities and the power production is mostly influenced by the water flow of the river. Reservoir power plants are located beside a dam, which is used to control the volume of water in the reservoirs to modulate power production. The solution to the short-term scheduling problem determines the optimal unit commitment and loading, the water flow through the turbines and the volume of water in the reservoirs, for every period, usually hourly or daily. This problem is complex to solve since the hydropower production functions are nonconvex and nonlinear. Also, each turbine has a different efficiency and it is important to limit unit restarts since they cause premature wear to the equipment. Finally, precipitations that fall on the watersheds are uncertain, which means that the reservoir inflows are unknown when a decision is taken. This thesis proposes a modeling of the short-term problem for the Saguenay Lac-St-Jean hydroelectric system, and proposes a stochastic optimization model to solve it.

First, a deterministic modeling of the problem is proposed, which means that the inflows are
known. A dynamic programming algorithm is used to calculate power production functions
for each power plant, which are then used in the objective functions of a two-phase opti-
mization process. The first phase, namely the loading problem, is a nonlinear mixed integer
program. The second phase, the unit commitment, is a linear integer model. Results obtained
from the two phase optimization process are compared to historical operational decisions and
show that the method increases energy production. Also, a comparison with a modeling of
the problem in a single optimization model confirms that the proposed methodology allows
to solve rapidly the short-term unit commitment and loading problem.

Second, the deterministic modeling is extended to consider uncertain water inflows in the reservoirs and a stochastic short-term model is developed. A scenario tree generation method, based on the minimization of the nested distance, is used to create scenario tree from the inflow forecasts that are then used as input to the stochastic optimization model. The deterministic equivalent of the model is solved, which means that a deterministic model is solved at every node of the scenario tree. A rolling-horizon simulation is used to validate the optimization models. Inflow forecasts are updated daily, a scenario tree is generated, the stochastic optimization model is solved, the solution of the first scenario tree node is implemented, reservoir volumes are updated with the real realizations of the inflows and the process is relaunched for the whole planning horizon. Numerical results compare the solution obtained from the stochastic optimization to the solution obtained with the use of the median scenario of the inflows forecast only and show that the stochastic model allows to produce more energy. Multiple numerical tests are conducted, with different numbers of stages and scenarios in the scenario trees.

Finally, the complexity required in the structure of the scenario trees to maximize energy pro-
duction in a rolling-horizon framework is investigated. The scenario tree generation method
requires input parameters which are the number of stages, the number of child nodes at each
stage and the aggregation of the period covered by each stage. The best scenario tree param-
eters are found using a Blackbox optimization formulation of the stochastic short-term unit
commitment and loading problem that maximizes the energy production over the rolling-
horizon. Three comparisons with the stochastic short-term model are conducted. The first
one involves generating a set of scenario trees built from inflow forecast data over a rolling-
horizon. The second replaces the set of scenario trees by the median scenario. The last one
replaces the set of trees by scenario fans. Numerical results show that using a set of scenario
trees is better than using the median scenario, but using scenario fans yields a comparable
solution to using scenario trees with less computational effort.
TABLE DES MATIÈRES

<table><tr><td>DÉDICACE</td><td>iii</td></tr><tr><td>REMERCIEMENTS</td><td>iv</td></tr><tr><td>RÉSUMÉ</td><td>vi</td></tr><tr><td>ABSTRACT</td><td>viii</td></tr><tr><td>TABLE DES MATIÈRES</td><td>x</td></tr><tr><td>LISTE DES TABLEAUX</td><td>xiii</td></tr><tr><td>LISTE DES FIGURES</td><td>xiv</td></tr><tr><td>LISTE DES SIGLES ET ABRÉVIATIONS</td><td>xvi</td></tr><tr><td>CHAPITRE 1 INTRODUCTION</td><td>1</td></tr><tr><td>1.1 Définitions et concepts de base</td><td>1</td></tr><tr><td>1.2 Éléments de la problématique</td><td>3</td></tr><tr><td>1.3 Objectifs de recherche</td><td>5</td></tr><tr><td>1.4 Plan de la thèse</td><td>7</td></tr><tr><td>CHAPITRE 2 REVUE DE LITTÉRATURE</td><td>8</td></tr><tr><td>2.1 Méthodes déterministes</td><td>8</td></tr><tr><td>2.1.1 Programmation dynamique</td><td>9</td></tr><tr><td>2.1.2 Programmation mixte en nombres entiers</td><td>10</td></tr><tr><td>2.1.3 Relaxation lagrangienne</td><td>12</td></tr><tr><td>2.1.4 Autres approches</td><td>13</td></tr><tr><td>2.2 Méthodes stochastiques</td><td>13</td></tr><tr><td>2.2.1 Programmation dynamique stochastique</td><td>13</td></tr><tr><td>2.2.2 Méthodes stochastiques pour résoudre le problème à court terme</td><td>15</td></tr><tr><td>2.3 Optimisation boîte noire</td><td>21</td></tr><tr><td>2.3.1 Algorithme de recherche par coordonnées (CS)</td><td>21</td></tr><tr><td>2.3.2 Algorithme de recherche par motifs (PS)</td><td>22</td></tr><tr><td>2.3.3 Algorithme de recherche sur treillis adaptatifs (MADS)</td><td>22</td></tr><tr><td>2.3.4 Applications des méthodes d'optimisation boîte noire</td><td>24</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>CHAPITRE 6 ARTICLE 3 : SCENARIO TREE MODELING FOR STOCHASTIC SHORT-TERM HYDROPOWER OPERATIONS PLANNING</td>
      <td>76</td>
    </tr>
    <tr>
      <td>6.1 Introduction</td>
      <td>76</td>
    </tr>
    <tr>
      <td>6.2 Methodology</td>
      <td>79</td>
    </tr>
    <tr>
      <td>6.2.1 Stochastic Short-Term Unit Commitment And Loading Problem</td>
      <td>79</td>
    </tr>
    <tr>
      <td>6.2.2 Rolling-Horizon Scheme</td>
      <td>79</td>
    </tr>
    <tr>
      <td>6.2.3 Parameters Of The Scenario Tree</td>
      <td>81</td>
    </tr>
    <tr>
      <td>6.2.4 Loading And Unit Commitment Problems</td>
      <td>82</td>
    </tr>
    <tr>
      <td>6.3 Blackbox optimization</td>
      <td>83</td>
    </tr>
    <tr>
      <td>6.3.1 Blackbox Optimization Concept</td>
      <td>83</td>
    </tr>
    <tr>
      <td>6.3.2 Blackbox Formulation Of The Rolling-Horizon Scheme</td>
      <td>83</td>
    </tr>
    <tr>
      <td>6.3.3 Blackbox Validation</td>
      <td>86</td>
    </tr>
    <tr>
      <td>6.4 Case study</td>
      <td>88</td>
    </tr>
    <tr>
      <td>6.4.1 Hydropower System</td>
      <td>88</td>
    </tr>
    <tr>
      <td>6.4.2 Structure Of The Data</td>
      <td>89</td>
    </tr>
    <tr>
      <td>6.5 Computational experiments</td>
      <td>89</td>
    </tr>
    <tr>
      <td>6.5.1 Deterministic Optimization With The Median Scenario</td>
      <td>90</td>
    </tr>
    <tr>
      <td>6.5.2 Scenario Fans</td>
      <td>92</td>
    </tr>
    <tr>
      <td>6.6 Conclusion</td>
      <td>94</td>
    </tr>
    <tr>
      <td>6.7 Notation</td>
      <td>95</td>
    </tr>
    <tr>
      <td>CHAPITRE 7 DISCUSSION GÉNÉRALE</td>
      <td>96</td>
    </tr>
    <tr>
      <td>CHAPITRE 8 CONCLUSION ET RECOMMANDATIONS</td>
      <td>98</td>
    </tr>
    <tr>
      <td>RÉFÉRENCES</td>
      <td>101</td>
    </tr>
  </tbody>
</table>
LISTE DES TABLEAUX

<table><tr><td>Table 4.2</td><td>Turbine combinations at Chute-du-Diable</td><td>34</td></tr><tr><td>Table 4.3</td><td>Total energy production with same initial solution</td><td>43</td></tr><tr><td>Table 4.4</td><td>Benchmarking on instances with 10 turbines</td><td>47</td></tr><tr><td>Table 5.1</td><td>Results for 6 test cases (5 are data sets from the year 2011 and 1 from 2010). Energy produced by the stochastic solution and the median scenario rolling-horizon is given. Also, the difference in energy between both solutions is shown.</td><td>69</td></tr><tr><td>Table 5.2</td><td>Objective function values for 6 random scenario trees with the same number of stages and scenarios, on two data sets.</td><td>70</td></tr><tr><td>Table 6.1</td><td>Reservoir and power plants characteristics</td><td>89</td></tr></table>
LISTE DES FIGURES

<table>
  <tbody>
    <tr>
      <td>Figure 1.1</td>
      <td>Vue en coupe d'une centrale hydroélectrique. Adaptée avec permission de "Hydroelectric dam" par Tomia, 2008. Image sous license GFDL et CC-BY-2.5.</td>
      <td>2</td>
    </tr>
    <tr>
      <td>Figure 1.2</td>
      <td>Efficacité en fonction du débit turbiné unitaire</td>
      <td>4</td>
    </tr>
    <tr>
      <td>Figure 1.3</td>
      <td>Réseau hydroélectrique du Saguenay Lac-St-Jean. Les rectangles représentent des centrales et les triangles des réservoirs.</td>
      <td>6</td>
    </tr>
    <tr>
      <td>Figure 2.1</td>
      <td>Arbre de scénarios avec 4 scénarios et 4 étapes.</td>
      <td>17</td>
    </tr>
    <tr>
      <td>Figure 2.2</td>
      <td>Algorithme de recherche par coordonnées</td>
      <td>22</td>
    </tr>
    <tr>
      <td>Figure 2.3</td>
      <td>Optimisation d'une boîte noire.</td>
      <td>23</td>
    </tr>
    <tr>
      <td>Figure 4.1</td>
      <td>Turbine efficiency as a function of water discharge.</td>
      <td>33</td>
    </tr>
    <tr>
      <td>Figure 4.2</td>
      <td>Power output at Chute-du-Diable.</td>
      <td>35</td>
    </tr>
    <tr>
      <td>Figure 4.3</td>
      <td>Power output at Chute-du-Diable.</td>
      <td>36</td>
    </tr>
    <tr>
      <td>Figure 4.4</td>
      <td>Power output with different number of active turbines.</td>
      <td>37</td>
    </tr>
    <tr>
      <td>Figure 4.5</td>
      <td>Histogram of differences of energy.</td>
      <td>44</td>
    </tr>
    <tr>
      <td>Figure 4.6</td>
      <td>Case 1 water discharge and water head.</td>
      <td>45</td>
    </tr>
    <tr>
      <td>Figure 4.7</td>
      <td>Case 6 water discharge and water head.</td>
      <td>46</td>
    </tr>
    <tr>
      <td>Figure 4.8</td>
      <td>Case 9 benchmark and two-phase approach water head.</td>
      <td>48</td>
    </tr>
    <tr>
      <td>Figure 5.1</td>
      <td>Building inflow scenarios from a 7 day deterministic precipitation forecast.</td>
      <td>55</td>
    </tr>
    <tr>
      <td>Figure 5.2</td>
      <td>Generation of a random path based on three available data paths of inflows. The generated value of inflow for stage 1 is shown with a star marker.</td>
      <td>56</td>
    </tr>
    <tr>
      <td>Figure 5.3</td>
      <td>Illustration of the 2 steps of the algorithm. Generation of a random path of inflows from available data paths of inflows and stochastic approximation to improve the value of some scenario tree nodes.</td>
      <td>59</td>
    </tr>
    <tr>
      <td>Figure 5.4</td>
      <td>A scenario tree with node probabilities (over the node) and scenario probabilities (indicated with π<sub>i</sub>)</td>
      <td>60</td>
    </tr>
    <tr>
      <td>Figure 5.5</td>
      <td>Maximum output surfaces</td>
      <td>62</td>
    </tr>
    <tr>
      <td>Figure 5.6</td>
      <td>Hydroelectric system studied.</td>
      <td>66</td>
    </tr>
    <tr>
      <td>Figure 5.7</td>
      <td>Average computational time of scenario tree generation and optimization for one day in the rolling-horizon.</td>
      <td>68</td>
    </tr>
    <tr>
      <td>Figure 5.8</td>
      <td>Water discharges and reservoir levels for the case June 2011, 5 stages, 16 scenarios.</td>
      <td>71</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>Figure 5.9</td>
      <td>Water discharges and reservoir levels for the case June 2011, 7 stages, 16 scenarios.</td>
      <td>72</td>
    </tr>
    <tr>
      <td>Figure 5.10</td>
      <td>Comparison of September (lower figures in each subfigure) and October (upper figures in each subfigure) day 1 data sets. The dashed lines are the minimum and maximum scenarios. The median scenario is the solid line. The actual realization of the inflows is the plus sign line.</td>
      <td>73</td>
    </tr>
    <tr>
      <td>Figure 6.1</td>
      <td>Rolling-horizon scheme. These steps are repeated for each day h of the rolling-horizon.</td>
      <td>80</td>
    </tr>
    <tr>
      <td>Figure 6.2</td>
      <td>A scenario tree with D = 4 stages</td>
      <td>82</td>
    </tr>
    <tr>
      <td>Figure 6.3</td>
      <td>The Blackbox optimization solution approach. Index h is the day in the rolling-horizon framework.</td>
      <td>85</td>
    </tr>
    <tr>
      <td>Figure 6.4</td>
      <td>Validation of the blackbox. For the same scenario tree structure, histogram of the energy production in GWh.</td>
      <td>86</td>
    </tr>
    <tr>
      <td>Figure 6.5</td>
      <td>Convergence of Blackbox optimization for a 5-stage scenario tree.</td>
      <td>87</td>
    </tr>
    <tr>
      <td>Figure 6.6</td>
      <td>Hydroelectric system studied. Squares represent power plants and triangles reservoirs.</td>
      <td>88</td>
    </tr>
    <tr>
      <td>Figure 6.7</td>
      <td>Scenario trees built from inflow forecast scenarios</td>
      <td>89</td>
    </tr>
    <tr>
      <td>Figure 6.8</td>
      <td>Boxplots of the difference in energy production between stochastic and deterministic optimization</td>
      <td>91</td>
    </tr>
    <tr>
      <td>Figure 6.9</td>
      <td>Scenario fan with x scenarios and D stages</td>
      <td>92</td>
    </tr>
    <tr>
      <td>Figure 6.10</td>
      <td>Boxplots of the difference in energy production between stochastic optimization of the structure and the aggregation of the scenario tree compared to scenario fans</td>
      <td>93</td>
    </tr>
  </tbody>
</table>
LISTE DES SIGLES ET ABRÉVIATIONS

ARIMA Auto-Regressive Integrated Moving Average

CS Coordinate Search

GPS Generalized Pattern Search

MADS Mesh Adaptive Direct Search

PS Pattern Search
# CHAPITRE 1 INTRODUCTION

La province de Québec est située sur un territoire riche en rivières, ce qui en fait un endroit de prédilection pour la production hydroélectrique. Cette province possède trois pour cent de la réserve mondiale en eau douce et dix pour cent de son territoire en est couvert. La production d'aluminium est un processus très énergivore, résultant du procédé de fusion qui permet de réduire l'alumine, obtenue à partir de bauxite, par électrolyse. Ce procédé consiste à faire fondre le mélange d'alumine et de cryolithe dans une cuve, puis à y faire traverser un courant électrique continu par une anode et une cathode. L'aluminium liquide se dépose ainsi au fond de la cuve où il peut être récupéré. Plusieurs producteurs d'aluminium sont installés au Québec à proximité de centrales hydroélectriques afin d'alimenter leurs usines en énergie hydroélectrique. C'est notamment le cas de Rio Tinto, une entreprise minière qui exploite des usines de production d'aluminium dans la région du Saguenay Lac-St-Jean. Cette compagnie est aussi propriétaire d'un système de production hydroélectrique lui permettant d'alimenter la demande en énergie de ses usines à hauteur de 90%. Ils doivent acheter l'énergie manquante à Hydro-Québec, donc il est primordial pour eux de gérer leur réseau hydroélectrique de façon efficace pour produire le plus d'énergie possible avec la quantité d'eau disponible afin de limiter leur achat d'énergie.

## 1.1 Définitions et concepts de base

La gestion d'un système de production hydroélectrique est complexe et nécessite plusieurs processus d'optimisation. Un tel système est constitué de centrales hydroélectriques qui abritent les groupes turbo-alternateurs, ou turbines, qui transforment l'énergie mécanique produite par la pression de l'eau sur les pales de la turbine en énergie électrique. Deux types de centrales existent : au fil de l'eau et à réservoir. Les centrales au fil de l'eau ont peu de réserve d'eau et la production d'énergie est surtout influencée par le débit du cours d'eau. Les centrales à réservoir sont situées à côté d'un barrage qui permet de moduler la production d'énergie en variant la hauteur de chute, qui est la différence entre l'élévation de l'eau en amont et en aval du barrage. Une conduite forcée, aussi connue sous le nom de conduite d'amenée, relie le barrage à la centrale qui, comme son nom l'indique, permet d'amener l'eau à la turbine.

La figure 1.1 montre une vue en coupe d'une centrale hydroélectrique. La centrale hydroélectrique B est composée d'un groupe turbo-alternateur dénoté par C et D. Le bief amont, qui est la hauteur d'eau maximale du barrage est désigné par A. La conduite d'amenée, en F,
amène l'eau retenue par le barrage à la turbine et H est la rivière en aval du barrage, où le bief aval est la hauteur d'eau maximale à la sortie.

<figure><img src="image_1.png" /><figcaption>Figure 1.1 Vue en coupe d'une centrale hydroélectrique. Adaptée avec permission de "Hydroelectric dam" par Tomia, 2008. Image sous licence GFDL et CC-BY-2.5.</figcaption></figure>

Plusieurs facteurs influencent la capacité de production d'énergie des centrales hydroélectriques et le plus important est sans contredit les précipitations qui tombent sur le bassin hydrographique. Aussi désigné sous le nom de bassin versant, c'est le territoire sur lequel tombent les précipitations sous forme de pluie et de neige, nommés apports naturels des réservoirs. Ces données sont incertaines puisque les précipitations ne sont pas connues avant de prendre une décision pour l'opération du système hydroélectrique. Une bonne prévision des phénomènes naturels est primordiale afin d'assurer la production d'énergie. Des modèles mathématiques sont utilisés pour tenter de prévoir le plus fidèlement possible les apports naturels des réservoirs. D'autres modèles d'optimisation sont nécessaires pour répartir la production de façon optimale entre toutes les centrales et tous les groupes turbo-alternateurs d'un système de production hydroélectrique. Plusieurs modèles permettent d'effectuer cette gestion : long terme, moyen terme et court terme.

L'optimisation à long terme est utilisée pour déterminer le potentiel de production futur, considérant l'incertitude des apports naturels des réservoirs et son horizon est de quelques années. Elle n'est pas utilisée sur une base opérationnelle régulière, mais plutôt lors de modifications majeures au réseau.

L'optimisation moyen terme est utilisée sur une base hebdomadaire afin d'estimer la quantité d'eau disponible à la production. Ce modèle permet de gérer les réservoirs des centrales et
prend en considération l'écoulement de l'eau entre les centrales, les niveaux des réservoirs, les apports naturels ainsi que leur incertitude et la demande en énergie, pour ne nommer que ceux-là.

L'optimisation à court terme est nécessaire afin de déterminer comment répartir l'eau disponible à la production, obtenue de l'optimisation moyen terme, entre les turbines et les centrales. Son horizon est court, habituellement des périodes allant de une à 24 heures sur un horizon d'une à quelques semaines.

Cette thèse traite uniquement de l'optimisation à court terme.

## 1.2 Éléments de la problématique

Les problèmes d'optimisation à court terme sont complexes à résoudre, car les fonctions de production hydroélectrique sont non convexes. Elles sont aussi non linéaires et dépendent de l'efficacité des turbines, de la hauteur de chute nette qui est une fonction non linéaire du débit turbiné et du volume, et, finalement, du débit turbiné unitaire.

La puissance pouvant être générée par une turbine [106] est donnée par l'équation :

$$P = n(Q) \times g \times Q \times h_n(v, Q), \quad (1.1)$$

où P est la puissance en kW, n(Q) est l'efficacité de la turbine, g est l'accélération gravitationnelle en m/s², Q est le débit turbiné en m³/s, v est le volume du réservoir en hm³ et h_n est la hauteur de chute nette en m. Pour une turbine donnée, la puissance est une fonction du débit turbiné, de la hauteur de chute nette et de l'efficacité. La hauteur de chute brute est la différence entre le bief amont et le bief aval. Lorsque l'eau circule dans la conduite d'amenée, la friction cause des pertes sous forme de chaleur, réduisant l'énergie. Ce phénomène crée des pertes de charges qui doivent être considérées dans le calcul de la puissance, avec l'utilisation de la hauteur de chute nette, qui est la hauteur de chute brute à laquelle les pertes de charge sont soustraites. La hauteur de chute nette est calculée par :

$$h_n(v, Q) = h_f(v) - h_t(Q) - \psi(Q), \quad (1.2)$$

où v est le volume du réservoir en hm³, Q est le débit total à la centrale en m³/s, h_f est une fonction non linéaire retournant l'élévation amont en m, h_t est une fonction non linéaire retournant l'élévation aval en m et ψ est une fonction non linéaire retournant les pertes de charge en m.

Une des particularités du système étudié est qu'il n'y a pas de représentation analytique des
fonctions de puissance. D'autres caractéristiques du problème compliquent sa résolution et sont expliquées ci-après.

**Efficacité** Chaque turbine possède sa propre courbe d'efficacité, ce qui signifie qu'elles ne produiront pas la même puissance pour un débit unitaire et une hauteur de chute donnés. À titre d'exemple pour une turbine, la figure 1.2 montre que le débit turbiné unitaire maximal est de 48 m³/s. Au-delà de cette valeur, il est nécessaire de déverser de l'eau, ce qui fait diminuer la puissance et donc l'efficacité. Le rendement maximal de la turbine est obtenu pour un débit de 41 m³/s. Dans ces conditions, la turbine est opérée dans sa zone optimale, autant pour maximiser la puissance que la durée de vie de l'équipement.

<figure><img src="image_2.png" /><figcaption>Figure 1.2 Efficacité en fonction du débit turbiné unitaire</figcaption></figure>

**Zones d'opération interdites** Les turbines ont des zones d'opération interdites. Sous certaines conditions d'opération, un vortex peut se former dans la turbine et créer des variations de pression qui peuvent causer des dommages aux composantes mécaniques. Elles ne doivent pas être opérées sous ces conditions. Les zones interdites sont souvent discontinues lorsqu'une turbine en a plusieurs.

**Démarrages de turbines** Les démarrages fréquents de turbines réduisent leur durée de vie. Un démarrage peut être considéré comme un certain nombre d'heures d'opération. Ainsi, chaque fois qu'une turbine est démarrée, sa durée de vie est raccourcie. Cette caractéristique
est modélisée en terme de coût de démarrage dans les modèles d'optimisation. La valeur du coût est calculée par Rio Tinto et est utilisée comme un paramètre dans le modèle.

Valeur de l'eau L'objectif du modèle d'optimisation à court terme est de maximiser l'énergie produite. Dans un contexte d'application réel, il est crucial d'empêcher l'optimisation de vider le réservoir à la fin de l'horizon de planification puisqu'une continuité est nécessaire d'une période à l'autre. Si l'eau à la fin de l'horizon de planification n'est pas valorisée, l'optimisation videra le réservoir afin de turbiner toute l'eau qui s'y trouve, peu importe qu'elle le soit de façon efficace ou non puisque de l'énergie sera produite. Dans un modèle d'optimisation déterministe, c'est-à-dire avec un seul scénario d'apports, il y a deux façons de s'assurer que le réservoir ne sera pas vidé. La première est de contraindre le volume final du réservoir et la seconde est d'avoir une fonction valeur finale de l'eau qui valorise l'eau à la fin de l'horizon dans l'objectif. Il est possible d'obtenir le volume final du réservoir par un modèle d'optimisation à moyen terme, par exemple. Cela signifie que les apports sont supposés parfaitement connus tout au long de l'horizon de planification et ainsi le volume final du réservoir peut être obtenu et utilisé comme paramètre dans le modèle d'optimisation.

Or, avec un modèle d'optimisation stochastique, dans lequel plusieurs scénarios d'apports sont considérés, imposer une cible de volume final devient problématique puisque les apports naturels sont inconnus et que chaque scénario a une probabilité de se réaliser. Cela dit, un modèle d'optimisation à moyen terme peut être résolu pour tous les scénarios, mais chaque scénario a une valeur de volume final différente. Ainsi, une cible finale sur le volume pourrait ne pas être réalisable pour un autre scénario si l'effet des conditions hydriques est négligé, ce qui implique que la cible finale est propre à un scénario et que la solution trouvée ne sera possiblement pas optimale. Une façon de contraindre l'optimisation à ne pas vider le réservoir est de définir une fonction valeur de l'eau à la fin de l'horizon qui est maximisée. Cette fonction est construite avec un modèle d'optimisation moyen terme à l'aide d'un algorithme de programmation dynamique stochastique et elle dépend du volume du réservoir. Pour les besoins de cette thèse, cette fonction est élaborée par Rio Tinto et utilisée dans l'objectif du modèle d'optimisation à court terme.

## 1.3 Objectifs de recherche

Cette thèse, rédigée par articles, a pour objectif de développer un modèle d'optimisation stochastique à court terme pour le réseau hydroélectrique du Saguenay Lac-St-Jean. Le modèle considère toutes les caractéristiques du problème expliquées précédemment en plus de tenir compte de l'incertitude liée aux apports naturels des réservoirs. Plus spécifiquement, un
modèle d'optimisation stochastique à court terme est développé et permet de déterminer les volumes des réservoirs, les débits turbinés et les turbines en marche à chaque pas de temps, pour chaque centrale.

Le réseau hydroélectrique est composé de cinq centrales hydroélectriques soit : Chute-Savane, Chute-du-Diable, Shipshaw/Chute-à-Caron, Chute-des-Passes et Isle-Maligne, illustrées à la figure 1.3. Le réseau est composé de 42 turbines et la capacité installée est de 3100 MW. Quatre réservoirs sont disponibles et trois d'entre eux ont une capacité de stockage supérieure à 2000 $hm^3$. Pour les besoins de la thèse, le modèle d'optimisation à court terme est validé sur un sous-ensemble de centrales, soit Chute-du-Diable et Chute-Savane pour le premier article et en ajoutant Isle-Maligne pour les deuxième et troisième articles.

<figure><img src="image_3.png" /><figcaption>Figure 1.3 Réseau hydroélectrique du Saguenay Lac-St-Jean. Les rectangles représentent des centrales et les triangles des réservoirs.</figcaption></figure>

Le premier objectif consiste à développer un modèle d'optimisation déterministe à court terme pour le réseau hydroélectrique du Saguenay Lac-St-Jean. Cet objectif permet de modéliser le problème à court terme et la solution obtenue par l'optimisation est comparée aux décisions opérationnelles historiques de Rio Tinto. Le deuxième objectif reprend cette modélisation mais la modifie pour formuler un modèle d'optimisation stochastique, c'est-à-dire que les
apports naturels sont incertains. Une méthode de génération d'arbres de scénarios, qui permettent de représenter l'incertitude des apports dans le modèle d'optimisation à court terme, est utilisée. Finalement, le dernier objectif consiste à déterminer la complexité nécessaire dans les arbres de scénarios afin d'obtenir une bonne solution au modèle d'optimisation à court terme. Pour ce faire, l'optimisation boîte noire est utilisée pour optimiser les paramètres de la méthode de génération d'arbres de scénarios.

## 1.4 Plan de la thèse

Le chapitre deux présente une revue de littérature des principales méthodes d'optimisation court terme déterministes et stochastiques, ainsi que de l'optimisation boîte noire. Le troisième chapitre détaille l'organisation de la thèse ainsi que la démarche scientifique liant les articles. Les chapitres quatre, cinq et six présentent les articles. Une discussion générale synthétisant le travail effectué est présentée dans le chapitre sept. Finalement, la conclusion est exposée au chapitre huit, discutant des limitations de la solution proposée ainsi que des perspectives de recherche futures.
# CHAPITRE 2 REVUE DE LITTÉRATURE

Les notions présentées dans ce chapitre sont nécessaires à la bonne compréhension des trois objectifs de recherche présentés dans cette thèse. Le premier objectif est de créer un modèle à court terme déterministe. Le second est d'utiliser une méthode de génération d'arbres de scénarios afin de développer un modèle d'optimisation stochastique qui considère l'incertitude liée aux apports naturels des réservoirs. Le troisième objectif étudie la complexité nécessaire dans les arbres de scénarios utilisés dans le modèle à court terme, en utilisant l'optimisation boîte noire pour optimiser les paramètres de génération des arbres de scénarios.

Ce chapitre est découpé en plusieurs parties. Tout d'abord, le problème de chargement et de répartition déterministe est abordé. Par la suite, les méthodes d'optimisation stochastiques sont présentées. Finalement, une revue du concept d'"optimisation boîte noire" est effectuée.

## 2.1 Méthodes déterministes

Le premier objectif de recherche est de développer un modèle déterministe. Le problème de chargement et de répartition optimal, ou à court terme, est complexe et plusieurs interdépendances existent entre les variables qui composent ce problème. Le volume et le débit doivent être optimisés simultanément, ce qui complique sa résolution. De plus, tel qu'expliqué à l'introduction, chaque turbine possède une efficacité différente, les turbines ont des zones d'opération interdites, et il est important de limiter les démarrages. Dans ces modèles, il n'y a pas d'incertitude et ils sont déterministes [40, 96]. Les fonctions objectif des problèmes à court terme dépendent de la modélisation mais peuvent être de : maximiser la production d'énergie, minimiser les coûts d'opération ou minimiser les pertes de production d'énergie, par exemple. Les contraintes que l'on retrouve le plus souvent sont : les contraintes hydriques et les bornes sur les volumes/niveaux des réservoirs et les débits. Les contraintes hydriques représentent simplement la loi de la conservation de l'eau d'une période à l'autre.

Les principales méthodes de résolution, soient la programmation dynamique, la programmation mixte en nombres entiers et la relaxation lagrangienne, sont expliquées ci-après. Une revue récente [99] expose les différentes approches pour résoudre ce problème. Il est important de spécifier qu'au Québec, les producteurs doivent absolument transiger avec Hydro-Québec, ce qui signifie que les prix de l'énergie n'ont pas à être considérés dans les modèles d'optimisation, ce qui n'est pas le cas dans le reste du monde.
2.1.1 Programmation dynamique

Une des méthodes populaires pour résoudre le problème de chargement et de répartition est la programmation dynamique [24]. Cette méthode consiste à décortiquer le problème en sous-problèmes afin d'obtenir la solution optimale de façon récursive.

La programmation dynamique comporte des étapes et des états. Chaque sous-problème est une étape alors qu’à chaque étape, un état du système est défini. La solution est construite par induction arrière, c’est-à-dire que la résolution du problème commence à la dernière étape et les sous-problèmes sont résolus en reculant d’une étape à chaque fois, jusqu’à ce que tous les sous-problèmes soient résolus. Un de problèmes majeurs de cette méthode est que, rapidement, le problème devient difficile à résoudre compte tenu du nombre de variables. En fait, il est nécessaire de discrétiser l’espace d’états. Par exemple, en optimisation hydroélectrique, il faut discrétiser le volume, le débit unitaire, le débit total, sans oublier qu’il faut considérer le nombre de turbines ainsi que les périodes. Évidemment, cette énumération dépend de la modélisation du problème. Un des avantages de la programmation dynamique est que les non-linéarités sont traitées facilement, puisque toutes les combinaisons de variables sont évaluées.

Afin d’illustrer la méthode, la programmation dynamique est appliquée au problème de chargement des groupes, qui consiste à déterminer le débit turbiné par chaque turbine, selon un débit total $Q$ et un volume donné $v$. Chaque sous-problème correspond à une turbine et, en termes de programmation dynamique, une étape. Les états sont $s^j \in \{1, 2, ..., r\}$ le débit restant à turbiner pour la turbine $j$ et les variables d’optimisation sont $q^j \in \{1, 2, \min\{\bar{q}^j, Q\}\}$, le débit turbiné par la turbine $j$, où $\bar{q}^j$ est le débit maximum de la turbine $j$. L’objectif est de maximiser la puissance totale. La solution est retrouvée de façon récursive en résolvant un sous-problème par turbine. Pour $j = n$, où $n$ est le nombre total de turbines, la puissance optimale est $f^{*j}(s^j) = P(s^j, v)$, où $v$ est le volume donné et $P$ est la puissance. Le débit turbiné optimal est $q^{*j} = s^j$. Pour les autres sous-problèmes, c’est-à-dire pour $j = n-1, n-2, ..., 1$, la puissance optimale est donnée par $f^{*j}(s^j) = \max_{q^j} P(s^j, v) + f^{*j+1}(s^j - q^j)$ et le débit optimal par $q^{*j} = s^j$ qui maximise $f^{*j}(s^j)$.

Les auteurs de [7] utilisent la programmation dynamique pour déterminer le nombre de turbines en marche sur une base horaire. Pour ce faire, les étapes sont les heures, les états sont le nombre de turbines en marche à chaque étape et les variables de décision sont le nombre de démarrages et d’arrêts de turbines à chaque étape. Un compromis entre l’efficacité de la puissance produite et des démarrages est minimisé dans la fonction objectif du problème de programmation dynamique. Les auteurs ne considèrent pas les contraintes hydriques, car l’élévation amont ne varie pas beaucoup sur des périodes horaires. Pour des systèmes hy-
droélectriques comportant de petits réservoirs, cette hypothèse ne peut être appliquée. Une autre approche [109] utilise des approximations successives de la programmation dynamique, afin de contourner le problème de dimensionnalité. Une fonction de puissance représentant les pertes est utilisée afin de considérer l'efficacité des turbines, mais les élévations amonts sont négligées, ce qui est important à considérer lorsque les réservoirs ont des capacités de stockage relativement petites. Les étapes correspondent aux turbines, et il y a trois vecteurs d'état : la puissance accumulée des turbines, le débit turbiné total et la puissance maximale accumulée. Le but du problème est de maximiser l'efficacité totale du système de production hydroélectrique. Un autre article [80] utilise la programmation dynamique pour résoudre le problème de chargement et de répartition optimal dans le but de maximiser les revenus. Dans ce cas-ci, l'élévation amont est considérée, en plus des arrêts et démarrages des turbines. Une méthode en deux phases, comprenant chacune un algorithme de programmation dynamique est définie. La première phase construit un espace d'états et la seconde résout le problème de répartition et de chargement optimal. Cependant, le système hydroélectrique à l'étude comporte un seul réservoir. La modélisation considère les apports de réservoir, les prix de l'énergie et les disponibilités de turbines, mais l'ajout d'autres réservoirs et centrales impliquera des temps de calcul beaucoup plus élevés, voire rédhibitoires.

### 2.1.2 Programmation mixte en nombres entiers

Plusieurs auteurs proposent des méthodes de résolution du problème de chargement et de répartition à l'aide de modèles mixtes en nombres entiers. Ils peuvent être non linéaires [66] ou linéaires [49]. les fonctions de production hydroélectrique sont, soit dépendantes de la hauteur de chute et du débit turbiné, soit linéarisées et dépendent seulement du débit, en utilisant, par exemple, des fonctions linéaires par morceaux [31, 93] ou par interpolation. Ces problèmes comportent des variables continues pour le volume et le débit turbiné, par exemple, et des variables binaires pour l'arrêt et le démarrage des turbines.

Un des avantages de la formulation non linéaire est que l'effet de la hauteur de chute est considéré. De plus, les solveurs non linéaires actuels sont puissants et permettent de résoudre ces problèmes. Une formulation non linéaire de la fonction objectif qui dépend du débit turbiné et du volume est proposée et résolue avec un problème mixte non linéaire en nombres entiers [30]. Cette formulation est intéressante, car elle considère, entres autres, l'effet de la hauteur de chute et la pénalisation des démarrages de turbines. Une comparaison avec un modèle linéaire en nombres entiers, dont la fonction de production hydroélectrique varie linéairement avec le débit seulement, et ce modèle non linéaire est effectuée et les résultats montrent que la formulation non linéaire mixte permet d'obtenir de bons résultats dans un
temps de calcul satisfaisant, pour un système hydroélectrique de sept réservoirs et centrales. Pour un système de plus grande taille, les temps de calcul seraient beaucoup plus élevés. Les mêmes auteurs proposent une formulation quadratique du problème [29]. La fonction objectif comporte un terme pour maximiser les profits et un second pour valoriser l'eau contenue dans les réservoirs à la fin de l'horizon de planification. La formulation quadratique, comparativement à une formulation linéaire, permet d'augmenter les profits, mais il n'y a pas de comparaison avec une formulation non linéaire, ce qui aurait été intéressant, puisque la comparaison aurait été plus juste. Une autre approche consiste à minimiser les pertes de puissance [59] et à résoudre la relaxation du modèle non linéaire mixte en nombres entiers. À l'aide d'une méthode de simulation, la solution obtenue est ensuite ajustée afin de respecter les contraintes hydriques. Cette méthode est intéressante, car la décomposition permet de résoudre des problèmes de grand taille, dans ce cas-ci 95 centrales et 447 turbines, mais démontre que pour résoudre des problèmes réels, les formulations non linéaires mixtes en nombres entiers doivent être couplées à des méthodes de simulation et de décomposition, sinon elles sont impossibles à résoudre.

La modélisation est très importante dans les modèles linéarisés, car la solution peut être grandement influencée par la précision de celle-ci. Habituellement, l'effet de la hauteur de chute est implicite dans les modèles linéaires. Une technique d'interpolation [65] est utilisée afin de linéariser les fonctions de production hydroélectrique et considère l'effet de la hauteur de chute. Un problème linéaire mixte en nombres entiers est résolu afin de minimiser les coûts d'opérations du plus gros système de production hydroélectrique au monde, situé en Chine. Toutefois, cette modélisation ajoute un grand nombre de contraintes dans le modèle d'optimisation. Des tests sont réalisés pour un réservoir et 32 turbines. Advenant le cas où il y a plus de réservoirs, cette modélisation deviendrait rapidement difficile à appliquer, compte tenu du nombre de variables et de contraintes. Dans [23], les fonctions de production hydroélectrique sont linéarisées et l'effet de la hauteur de chute est considéré. Un modèle linéaire mixte en nombres entiers est résolu pour un système comportant un réservoir avec capacité de pompage. Un gros effort de linéarisation est réalisé et permet une bonne représentation des fonctions de production, mais tel qu'indiqué par les auteurs, cette technique est difficilement applicable en pratique, compte tenu de la taille du problème à résoudre.

Ces méthodes fonctionnent bien, mais plus le système hydroélectrique est complexe, plus elles deviennent lentes en termes de temps de calcul puisque le nombre de variables d'optimisation explose. Cependant, les solveurs modernes permettent de résoudre des problèmes de grande taille, mais plus le modèle est précis, plus difficiles ils sont à résoudre.
### 2.1.3 Relaxation lagrangienne

Les problèmes d'optimisation hydroélectrique à court terme deviennent rapidement difficiles à résoudre, car ils comportent un grand nombre de variables d'optimisation. De plus, selon la modélisation désirée, plusieurs inter-dépendances peuvent être considérées et compliquent la résolution des problèmes telles l'efficacité des turbines, l'effet de la hauteur de chute, les zones d'opération interdites ou les démarrages de turbines pour ne nommer que celles-là.

La relaxation lagrangienne [51, 98] permet de simplifier la résolution en séparant le problème original en plusieurs sous-problèmes faciles à résoudre. Les contraintes liantes, difficiles à traiter, sont relaxées et un problème maître coordonne la résolution en se basant sur la théorie de la dualité. Pour ce faire, les contraintes liantes, soit les contraintes hydriques, sont incorporées dans la fonction objectif avec un multiplicateur de Lagrange associé à chacune. Par la suite, le dual du problème lagrangien est résolu afin de trouver les valeurs des multiplicateurs. Cette façon de faire permet de trouver une borne inférieure serrée. Des méthodes de sous-gradients [76], qui sont des méthodes d'optimisation itératives, sont fréquemment utilisées pour résoudre le dual lagrangien.

Dans l'article [42], une méthode de décomposition en deux phases est proposée. La première phase est une relaxation lagrangienne, qui permet d'obtenir une borne sur la solution et la deuxième phase utilise une méthode de lagrangien augmenté [26] afin d'obtenir une solution réalisable. La méthode est validée sur un système hydroélectrique comportant 4 centrales et considère l'effet de la hauteur de chute, l'efficacité des turbines, les zones interdites et les démarrages des turbines sont pénalisés. Les résultats obtenus démontrent que cette méthode permet d'obtenir des résultats satisfaisants et que les temps de calcul sont raisonnables. Récemment, ces mêmes auteurs [43] ont comparé des méthodes de construction du dual lagrangien, appliquées au problème de répartition et de chargement optimal. Dans cet article, deux méthodes sont comparées : la relaxation des contraintes hydriques et l'ajout de variables artificielles qui sont relaxées. Afin de s'assurer que les solutions trouvées par le dual sont réalisable pour le primal, une méthode de lagrangien augmenté est utilisée, ce qui implique l'ajout d'un terme pour pénaliser les violations de contraintes dans la fonction objectif. Les turbines sont considérées identiques, ce qui n'est pas le cas en pratique. Les résultats montrent que la relaxation des contraintes hydriques permet d'obtenir de meilleurs résultats, mais qu'ils sont dépendants du problème. Ainsi, si le problème primal était différent, la conclusion pourrait aussi l'être, il est donc difficile de définir quelle méthode est la meilleure.

Une étude récente [44] compare les solutions obtenues par relaxation lagrangienne, programmation linéaire mixte en nombres entiers et programmation non linéaire mixte en nombres entiers. D'un point de vue opérationnel, la relaxation lagrangienne est la méthode qui requiert
le plus faible temps de calcul et qui peut résoudre des problèmes de grande dimension. Or,
son implémentation en pratique requiert beaucoup de connaissances et de programmation,
alors que les modèles mixtes en nombres entiers peuvent être utilisées plus facilement à l'aide
de solveurs commerciaux.

2.1.4 Autres approches

D'autres approches ont été utilisées afin de résoudre le problème de répartition et de char-
gement optimal. Par exemple, les algorithmes génétiques [79], des heuristiques basés sur les
fourmilières [70] ou des réseaux [78], mais ces méthodes requièrent une très bonne connais-
sance du problème pour l'ajustement de paramètres, ce qui n'est pas le cas lorsque des
méthodes mathématiques avec preuves de convergence, telles que présentées ci-haut, sont
utilisées. Toutefois, ces méthodes sont utiles pour traiter efficacement des problèmes non
convexes, non lisses ou discontinus.

2.2 Méthodes stochastiques

Le second objectif de recherche est de développer un modèle d'optimisation stochastique
à court terme pour le réseau hydroélectrique du Saguenay Lac-St-Jean. Cette section est
constituée de trois parties. La première expose la programmation dynamique stochastique, qui
est souvent utilisée pour élaborer la fonction de valorisation de l'eau restant dans le réservoir
à la fin de l'horizon de planification. La seconde partie présente les méthodes d'optimisation
stochastiques utilisées pour résoudre le problème d'optimisation à court terme. Finalement,
des méthodes de génération d'arbres de scénarios sont présentées.

2.2.1 Programmation dynamique stochastique

Lorsque des apports naturels incertains sont considérés dans le modèle d'optimisation à court terme, il est important de valoriser l'eau restant dans les réservoirs à la fin de l'horizon de planification afin d'empêcher les réservoirs de se vider. Comparativement à un modèle déterministe, un modèle stochastique ne peut utiliser de cible finale pour les volumes de réservoirs puisque le scénario d'apport qui se produira est inconnu. Il est possible de calculer le volume final pour chaque scénario, mais il est alors problématique de donner une seule valeur au modèle, puisque cette dernière pourrait ne pas être réalisable pour un scénario donné. L'optimisation moyen terme permet de gérer les réservoirs d'un système de production hydroélectrique. Elle est aussi utilisée en optimisation court terme déterministe afin de fournir le volume final du réservoir ou en optimisation court terme stochastique pour fournir une
fonction de valeur de l'eau au modèle à court terme. La programmation dynamique stochastique [18], qui est une extension de la programmation dynamique, est largement utilisée en optimisation moyen terme. D'autres méthodes découlent de celle-ci, par exemple la programmation dynamique stochastique duale [92] et la programmation dynamique stochastique par échantillonnage [32], mais pour les besoins de cette revue de littérature, seules les notions théoriques de programmation dynamique stochastique sont présentées. Tout comme pour la programmation dynamique, l'inconvénient majeur de la programmation dynamique stochastique est que les problèmes deviennent rapidement impossibles à résoudre, compte tenu de leur taille, surtout que l'incertitude doit maintenant être discrétisée.

La méthode de programmation dynamique stochastique appliquée à l'optimisation hydroélectrique moyen terme comporte plusieurs éléments :

- **Les états** à chaque étape sont notés $v_k^c$: le volume à la période $k$ et à la centrale $c$. Notons $S_k^c$ l'ensemble des états.

- **Les décisions** à chaque étape sont $q_k^c$: le débit turbiné à la période $k$ et à la centrale $c$. Notons $U_k^c(v_k^c)$ l'ensemble des décisions possibles selon l'état $v_k^c$.

- **Les variables aléatoires** sont $\delta_k^c$: l'apport naturel à la période $k$ pour la centrale $c$.

À chaque étape, selon l'état du système et la réalisation de la variable aléatoire, une puissance $\chi_k^c(v_k^c, q_k^c, \delta_k^c)$ est produite et elle est cumulative sur toutes les périodes. Il est important de noter que les apports des réservoirs sont corrélés dans le temps. Devant cette situation, l'introduction d'une variable hydrologique est nécessaire, afin de considérer les apports de la période précédente ou toute autre information permettant de bien représenter le système hydroélectrique à l'étude. Ainsi, la variable hydrologique $\psi_k^c$ peut être définie ainsi :

$$ \psi_k^c = \delta_{k-1}^c. \qquad (2.1) $$

Des politiques sont dites admissibles si une séquence de fonctions $\pi = \{\mu_0, \dots, \mu_{N-1}\}$ est telle que $q_k^c = q_k^c(v_k^c)$ et $q_k^c(v_k^c) \in U_k^c(v_k^c), \forall v_k^c \in S_k^c$, où $N$ est le nombre de périodes.

L'équation de Bellman, qui permet de calculer la valeur espérée de la puissance pour une politique admissible, dépend de l'état initial du système. Cependant, le principe d'optimalité de Bellman stipule que peu importe l'état initial, si $\pi^* = (\mu_0^*, \mu_1^*, \dots, \mu_{N-1}^*)$ est une politique optimale pour le problème complet, elle l'est aussi pour une politique tronquée $\pi_t^* = (\mu_t^*, \mu_{t+1}^*, \dots, \mu_{N-1}^*)$.

La fonction de récurrence, ou de Bellman, est alors $P_k^c(v_k^c, \psi_k^c)$ et correspond à la valeur espérée optimale de la puissance à l'étape $k$ à la fin si le système est à l'état $v_k^c$ à l'étape $k$ et
que la variable hydrologique est $\psi_k^c$. Les équations de Bellman à résoudre sont alors :

$$P_k^{c*}(v_k^c, \psi_k^c) = \mathbb{E}_{\delta_k^c | \psi_k^c} \left[ \max_{q_k^c \in U_k^c(q_k^c)} \left( \chi_k^c(v_k^c, q_k^c, \delta_k^c) + P_{k+1}^c(v_{k+1}, \psi_{k+1}^c) \right) \right], k = 0, 1, \dots, N-1. (2.2)$$

Dans [37], un modèle de programmation dynamique stochastique par échantillonnage est utilisé pour résoudre le problème d'optimisation des réservoirs. Comparativement à la programmation dynamique stochastique, cette méthode utilise des scénarios pour les apports naturels des réservoirs et les probabilités de transition pour passer d'un scénario i au scénario j doivent être calculées. Cela permet de préserver la persistance des apports d'une période à l'autre. La particularité de la méthode présentée dans l'article est l'utilisation des prévisions d'ensembles au lieu de l'historique des apports. L'avantage de ces prévisions, comparative aux données historiques, est qu'elles sont cohérentes avec les conditions actuelles du réservoir. De plus, la corrélation spatiale des apports entre plusieurs réservoirs est conservée puisque différentes réalisations tirées des bases de données historiques sont utilisées.

Tel que décrit dans [32], plusieurs choix de variable hydrologique sont possibles. Par exemple, l'apport de la période avant ou le meilleur apport prévu. L'introduction de cette variable dans le modèle permet de préserver les caractéristiques du système étudié, et le modèle donnera une approximation plus près de la réalité. En Amérique du Nord, l'indicateur de neige peut être un bon choix de variable hydrologique, puisque c'est cette variable qui permet d'approximer l'apport d'eau dans les bassins lors de la crue du printemps. Dans l'Ouest canadien et au Saguenay Lac-St-Jean, cette variable prend tout son sens puisque ce sont des régions qui reçoivent beaucoup de précipitations sous forme de neige.

Le défi de la programmation dynamique stochastique n'est pas d'utiliser l'algorithme en soi, qui est une extension de la programmation dynamique déterministe, mais plutôt de choisir la façon dont sera traitée l'incertitude. Selon la taille du système hydroélectrique à l'étude, il se peut qu'une modélisation déterministe utilisant la moyenne des apports naturels sur quelques années soit suffisante afin de résoudre le problème, tel que présenté dans [95], mais c'est rarement le cas dans des systèmes réels, où une bonne connaissance de l'historique permet de prendre de meilleures décisions.

### 2.2.2 Méthodes stochastiques pour résoudre le problème à court terme

Peu d'articles traitent du problème stochastique à court terme avec apports incertains dans la littérature. En considérant l'incertitude des apports, les turbines sont utilisées de façon plus efficace puisque l'optimisation stochastique est un compromis entre de faibles et de fortes prévisions d'apports. Par exemple, si les réservoirs sont pleins, le fait de considérer des
apports incertains lorsque des apports élevés sont prévus empêche l'optimisation de diminuer
le volume d'eau contenu dans les réservoirs et de forcer les turbines dans des zones non
efficaces, ce qui se traduit par des pertes de production d'énergie dans le futur si ces apports
élevés ne se produisent pas.

La majorité des articles qui traitent de l'optimisation court terme stochastique [3] utilisent
des arbres de scénarios pour représenter la distribution des apports incertains. Des modèles
de programmation stochastiques, résolus selon leur équivalent déterministe, c'est-à-dire que
chaque noeud de l'arbre de scénarios est un modèle d'optimisation déterministe, sont en-
suite résolus. Cette section présente les modèles de programmation stochastique ainsi que
différentes méthodes pour générer des arbres de scénarios.

Programmation stochastique avec recours

Un modèle de programmation stochastique comporte des paramètres incertains ou aléatoires.
La solution à ce problème doit donc être optimale pour toutes les valeurs que peuvent prendre
ce paramètre. Dans le modèle à court terme qui nous intéresse, les apports naturels sont
incertains. Les programmes stochastique peuvent être en deux étapes ou multi-étapes. Les
problèmes en deux étapes impliquent deux étapes de décisions. La première étape consiste à
prendre une décision avant de connaître la vraie réalisation de l'apport, alors que la deuxième
étape consiste à ajuster la décision prise à l'étape précédente lorsque la vraie réalisation de
l'apport est connue. Les problèmes multi-étapes sont simplement une extension des modèles
à deux étapes dans lesquels on considère des apports à plusieurs moments : lorsqu'un apport
est révélé, la décision correspondante est ajustée, une nouvelle décision est prise à l'étape
suivante et ainsi de suite. La référence [19] explique la programmation stochastique en détails,
veuillez vous y référer pour plus d'information. En pratique, les équivalents déterministes sont
résolus, ce qui signifie qu'un modèle d'optimisation est résolu pour toutes les possibilités, ou
les scénarios, possibles d'apports. Un arbre de scénarios est une représentation discrète d'un
processus stochastique continu et peut être vu à la figure 2.1. Cet arbre de scénarios a 4
scénarios et les probabilités que chacun se réalisent sont dénotées par π. Les probabilités de
chaque noeud constituant un scénario sont multipliées pour obtenir la probabilité du scénario.
Par exemple, le premier scénario est constitué des noeuds 1, 2, 3 et 5 et a une probabilité de
0,15 de se réaliser.
<figure><img src="image_4.png" /><figcaption>Figure 2.1 Arbre de scénarios avec 4 scénarios et 4 étapes.</figcaption></figure>

Cet arbre est ensuite utilisé comme paramètre d'entrée dans le problème d'optimisation stochastique multi-étape suivant, comportant une seule centrale hydroélectrique :

$$ \max_{q,v} \chi(v_1, q_1, \delta_1) + \sum_{j \in K} \pi_j \left( \sum_{i \in N_j} \chi(v_i, q_i, \delta_i) \right) + \sum_{t \in K} \pi_t \left( \sum_{p \in E_t} \Phi_p(v_p) \right) \quad (2.3) $$

$$ \text{sujet à :} \quad \delta_i = v_{i+1} - v_i + q_i, \quad \forall i \in N_j, \forall j \in K, \quad (2.4) $$

$$ v_{min} \le v_i \le v_{max}, \quad \forall i \in N, \quad (2.5) $$

$$ q_{min} \le q_i \le q_{max}, \quad \forall i \in N, \quad (2.6) $$

$$ v_i, q_i \ge 0, \quad \forall i \in N, \quad (2.7) $$

où K est l'ensemble des scénarios, N est l'ensemble des noeuds, E est l'ensemble des noeuds feuilles, $N_j$ est le sous-ensemble des noeuds du scénario $j$, $E_t$ est le sous-ensemble du noeud feuille pour le scénario $t$, $v_i$ est le volume, $q_i$ est le débit turbiné, $\delta_i$ est l'apport pour le noeud $i$ du scénario $N_j$, $\chi(v_1, q_1, \delta_1)$ est une fonction retournant l'énergie produite au premier noeud, $\Phi_p(v_p)$ est la fonction valeur de l'eau pour le noeud $p \in E_t$ et $\pi_j$ est la probabilité du scénario $j$. Noter que l'indice $i+1$ de la variable $v_{i+1}$ prend la valeur du noeud dans l'ensemble $N_j$. Par exemple, si $N_j = \{1, 3, 5, 7\}$ et que $i = 1, i+1 = 3$, car l'indice prend la valeur du noeud à la position $i = 2$ de l'ensemble $N_j$.

Le débit et le volume doivent être optimisés sans connaître le scénario d'apport qui se réalisera. La fonction objectif (2.3) consiste à maximiser l'énergie produite à la première étape et
l'espérance de la production d'énergie future, qui consiste en la somme de production d'énergie de tous les scénarios, multiplié par leur probabilité de réalisation. Aussi, le terme de droite permet de valoriser l'eau restante dans les réservoirs à la fin de l'horizon de planification. Les contraintes (2.4) sont celles de la conservation de l'eau et les volumes/débits minimums et maximums sont représentés par les contraintes (2.5) et (2.6). Finalement, la non-négativité est assurée par les contraintes (2.7).

Les programmes stochastiques sont non anticipatifs, ce qui signifie que les décisions prises à une certaine étape ne doivent pas être influencées par le futur, mais plutôt par les réalisations passées. La méthode de résolution habituelle est la création d'arbres de scénarios [35] afin d'approximer les apports suivant une distribution continue en les discrétisant.

Un des rares articles trouvé sur l'optimisation à court terme stochastique [45] considère les coûts de l'énergie et les apports naturels comme des variables incertaines. Une méthode basée sur celle-ci pourrait être envisagée, mais des modifications au modèle devraient être réalisées puisque le problème est situé dans un contexte de marché de l'énergie. L'incertitude est approximée par un arbre de scénarios de sept étapes et chaque étape est une journée divisée en 24 heures. L'objectif permet de trouver un équilibre entre le plan de production d'une journée et les profits futurs espérés. Dans [104], dix scénarios sont construits pour simuler l'incertitude sur les prix du marché et dix autres pour les apports naturels des réservoirs. Les données historiques et un modèle auto-régressif ARIMA (Auto-Regressive Integrated Moving Average) sont utilisés pour construire les prévisions. Dans le cadre de notre projet, les prévisions d'ensemble des apports naturels sont disponibles. L'étude montre que le modèle stochastique permet au producteur d'augmenter ses profits et que l'incertitude sur les prix permet un plus grand profit que l'incertitude sur les apports naturels. Tel qu'expliqué dans [27], les modèles peuvent être résolus de façon déterministe afin de diminuer l'effort de calcul, mais ils dépendent beaucoup du scénario. Ainsi, une solution optimale pour un scénario ne le sera pas nécessairement pour un autre. Un arbre de scénarios est alors utilisé afin de représenter l'incertitude et une probabilité de réalisation identique est donnée à chaque scénario. Une relaxation lagrangienne est effectuée sur les contraintes de non-anticipativité et sont pénalisées par la fonction objectif lorsqu'elles sont violées. La méthode donne de bons résultats, mais le choix du paramètre de pénalité est très important puisqu'il modifie la qualité de la solution, en plus du nombre d'itérations de l'algorithme. Les auteurs de [83] considèrent que la demande est incertaine et utilisent des arbres de scénarios pour considérer différents scénarios de demande. L'effet de la hauteur de chute est négligé, ce qui peut poser problème pour des systèmes qui comportent plusieurs réservoirs de faible capacité, même avec de courtes périodes.
## Génération d'arbres de scénarios

La génération des arbres de scénarios utilisés comme paramètres d'entrée au modèle d'optimisation stochastique est primordiale afin d'obtenir une bonne solution. En pratique, les données disponibles pour les apports sont souvent énormes et il est nécessaire de générer un arbre de scénarios basé sur ces données afin d'obtenir une représentation correcte de la distribution des apports et de pouvoir l'utiliser dans un modèle d'optimisation.

Plusieurs méthodes existent pour obtenir un arbre de scénarios à partir de données d'apports. Une revue des méthodes permettant de générer des scénarios est présentée dans l'article [61]. Les auteurs de [34] utilisent un modèle auto-régressif pour s'ajuster sur les données historiques des prix de l'énergie afin de générer des prix pour le modèle stochastique. L'arbre de scénarios est construit en échantillonnant la distribution du modèle ajusté. Une des difficultés d'utilisation de ce modèle auto-régressif est le choix des paramètres, car ils ont une grande influence sur les prix générés subséquemment. De plus, ils sont sensibles aux valeurs aberrantes, donc il est important de traiter les données d'entrée avant de bâtir le modèle.

D'autres méthodes créent une représentation discrète de la distribution des apports afin qu'elle corresponde à des propriétés statistiques pré-définies. Dans l'article [56], les quatre premiers moments, qui sont la moyenne, la variance, l'asymétrie et l'aplatissement de la représentation discrète correspondent aux quatre moments des scénarios d'apports initiaux de la distribution continue. Toutefois, cette méthode n'a pas de preuve de convergence et le nombre de scénarios nécessaires pour atteindre cette correspondance n'a pas de valeur fixe, ce qui fait en sorte qu'il est difficile d'évaluer le nombre de scénarios nécessaires, même si les quatre premiers moments sont identiques. Une autre méthode [82] minimise la distance imbriquée, qui est une généralisation de la distance de Wasserstein. Cette méthode itère deux phases jusqu'à ce que la convergence de la distance soit atteinte. La première phase échantillonne à partir de tous les scénarios d'apports initiaux et la deuxième phase utilise la méthode du gradient stochastique afin d'améliorer les valeurs d'apports de certains noeuds de l'arbre de scénarios. Cette méthode offre l'avantage d'utiliser tous les scénarios d'apports initiaux afin d'améliorer les valeurs dans l'arbre de scénarios et les quatre premiers moments correspondent toujours à la distribution continue des apports initiaux.

Une autre revue [35] présente plusieurs approches pour générer un arbre de scénarios incluant la recombinaison de données, les méthodes de contamination et aussi l'ajustement à des propriétés statistiques. Récemment, la méthode des copulas [60] a été utilisée pour générer des arbres de scénarios pour des arbres à deux étapes. Cette méthode offre l'avantage de traiter les dépendances mieux qu'en considérant seulement leur corrélation.
La méthode de réduction de scénarios [38, 52] est largement utilisée. Un arbre de scénarios initial est requis et une sélection directe, ou une réduction inverse est appliquée afin de réduire le nombre de scénarios et ainsi diminuer le temps de résolution du modèle d'optimisation stochastique. L'effet de la réduction sur la qualité de la solution, appliqué à un système de réservoirs hydroélectriques en cascade est disponible en [108].

## Autres approches

D'autres approches pour résoudre le modèle d'optimisation à court terme stochastique sont disponibles. L'optimisation robuste [13] et la programmation par contraintes probabilistes [103], par exemple. L'optimisation robuste résout des modèles stochastiques qui ont des paramètres incertains dans des ensembles incertains définis. Ainsi, l'optimisation cherche à trouver une solution qui est réalisable peu importe la réalisation du paramètre incertain. Dans [6] un horizon roulant est utilisé et l'optimisation robuste est appliquée aux décisions de la première journée alors que le reste de l'horizon est considéré déterministe. Cette méthode est intéressante puisque l'incertitude est appliquée aux décisions importantes. Un des problèmes de l'optimisation robuste est de définir l'incertitude. Dans les données historiques accessibles aux producteurs, certains scénarios d'apports peuvent être très faibles et d'autres très élevés. Ainsi, il est difficile de déterminer quelles sont les meilleures bornes pour les ensembles d'incertitude, en plus de capturer les non-linéarités présentes. Dans la programmation par contraintes probabilistes, les contraintes doivent être respectées selon une certaine probabilité. Un exemple appliqué à un réseau hydroélectrique en cascades se trouve dans [103]. Comme pour l'optimisation robuste, cette méthode exige des paramètres et probabilités en entrée, ce qui est difficile à faire en pratique.

## Horizon roulant

Dans un contexte opérationnel, les opérateurs des centrales hydroélectriques reçoivent des prévisions d'apports quotidiennement. Dès qu'elles sont disponibles, ils doivent lancer un modèle d'optimisation court terme afin de prendre une décision opérationnelle. Advenant qu'ils utilisent un modèle stochastique avec un arbre de scénarios pour représenter l'incertitude, seule la décision du premier noeud dans l'arbre de scénarios est utilisée. Des bancs d'essai en horizon roulant [74] sont utilisés pour imiter cette façon de faire pour valider les solutions des modèles d'optimisation stochastiques.

La méthodologie est la suivante. Puisque les prévisions d'apports sont mises à jour quotidiennement, un arbre de scénarios est généré, puis les modèles d'optimisation sont résolus et seule la solution du premier noeud de l'arbre de scénarios est implémentée. Lorsque la
vraie réalisation de l'apport est révélée, à la fin de la journée, les volumes des réservoirs sont mis à jour avec cette réalisation et la génération d'arbre de scénarios, ainsi que l'optimisation sont répétées pour la journée suivante, avec les apports mis à jour. Dans certains cas [17], la fenêtre de l'horizon de planification diminue lorsque les apports sont révélés et dans d'autres la fenêtre de l'horizon se déplace dans le temps [110]. Les auteurs de [33] utilisent un horizon roulant pour comparer quatre méthodes d'optimisation stochastiques. Les scénarios d'apports sont mis à jour à chaque fois qu'une décision est prise. La même méthodologie est retenue dans cette thèse.

## 2.3 Optimisation boîte noire

Le troisième objectif de recherche consiste à étudier la sensibilité de la structure des arbres de scénarios sur la gestion des centrales hydroélectriques dans un banc d'essai en horizon roulant. L'optimisation boîte noire est utilisée pour identifier les valeurs des paramètres de la méthode de génération d'arbres de scénarios qui permet d'obtenir une structure d'arbre qui maximise l'énergie produite dans l'horizon roulant. L'optimisation boîte noire permet de traiter des problèmes pour lesquels la fonction objectif ou les contraintes peuvent seulement être calculés par un code informatique. Les problèmes de boîtes noires sont souvent non lisses, non convexes et discontinus.

Le problème non contraint le plus simple est donné par :

$$ \min_{x \in \mathbb{R}^n} f(x), \qquad (2.8) $$

où f est une boîte noire.

L'algorithme MADS [10] (Mesh Adaptvie Direct Search) a été développé pour traiter des problèmes de boîtes noires. Deux principales classes d'algorithmes existent pour traiter ces problèmes : les méthodes de recherche directe, telle MADS, et les méthodes basées sur les modèles [84]. Dans cette thèse, nous présentons seulement les méthodes de recherche directe puisque ce sont celles qui seront utilisées. Pour les besoins de cette thèse, l'algorithme de recherche par coordonnées (CS), l'algorithme PS (Pattern Search) et MADS sont présentés.

### 2.3.1 Algorithme de recherche par coordonnées (CS)

L'ancêtre des méthodes de recherche directe, datant de 1952, est l'algorithme de recherche par coordonnées [39]. La figure 2.2 illustre la méthode. L'espace d'états est discrétisé en un treillis. La fonction objectif est évaluée pour un point initial, tel qu'illustré à la figure 2.2a. Ensuite,
les points de coordonnées nord, sud, est et ouest à ce point initial sont évalués, comme dans la figure 2.2b. Si un point permet d'abaisser la valeur de la fonction objectif, il est conservé et défini comme meilleure solution, tel que montré à la figure 2.2c. Advenant le cas où aucun point de coordonnée ne permet de diminuer la valeur de la fonction objectif, le treillis est réduit et la méthode est relancée, tel qu'illustré à la figure 2.2d. Lorsqu'une itération est un succès, c'est-à-dire qu'un nouveau point diminue la valeur de la fonction objectif, la taille du treillis est augmentée. L'algorithme s'arrête selon un critère d'arrêt pré-défini.

<figure><img src="image_5.png" /><figcaption>(a) Point initial</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) Points de coordonnées</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) Le point $f^*$ améliore la fonction objectif et est retenu</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(d) Le treillis est réduit</figcaption></figure>

<caption>Figure 2.2 Algorithme de recherche par coordonnées</caption>

## 2.3.2 Algorithme de recherche par motifs (PS)

L'algorithme PS [55] utilise deux types de déplacements. Le premier se nomme *exploratoire*. À partir du point initial, une coordonnée à la fois est modifiée et si la valeur de la fonction objectif est améliorée, la modification est conservée. Le deuxième type de déplacement est le *motif*. Comme son nom l'indique, les déplacements conservés dans la phase *exploratoire* sont appliqués à la solution courante, et le motif est appliqué de nouveau à la nouvelle solution. Si la valeur de la fonction objectif est améliorée avec cette nouvelle solution, alors le motif est répété sur cette solution. Sinon, l'algorithme retourne à la phase *exploratoire* pour trouver un nouveau motif de déplacement. Une extension de cet algorithme, nommée GPS (Generalized Pattern Search) [102] propose des modifications à la phase *exploratoire*.

## 2.3.3 Algorithme de recherche sur treillis adaptatifs (MADS)

L'algorithme MADS [10] est un algorithme qui permet de traiter des problèmes d'optimisation boîte noire. MADS est une extension des méthodes de recherche directes présentées ci-haut. Le domaine réalisable est discrétisé en un treillis dont la taille est ajustée à la fin de chaque itération. L'algorithme consiste en deux étapes qui sont répétées jusqu'à un critère d'arrêt pré-défini. La première étape, soit la phase de recherche, évalue différents points situés
sur le treillis afin de trouver une solution meilleure que la solution courante. Cette étape est flexible et peut être ajustée selon les spécificités du problème. La deuxième phase, soit l'étape de sonde, est obligatoire puisque les propriétés de convergence de MADS en découlent. Un ensemble générateur positif de directions est déterminé et si une solution meilleure que la solution courante est trouvée, elle est définie comme solution courante. À la fin de chaque itération de l'algorithme, le treillis est réduit lorsque la solution n'est pas améliorée et augmenté lorsqu'une nouvelle solution est trouvée. Cet algorithme est implémenté dans le logiciel NOMAD [63]. La différence majeure entre MADS et l'algorithme GPS est que MADS a deux paramètres pour le treillis, un pour la taille du treillis et un pour la taille de la sonde, comparativement à un seul pour GPS. Cette différence assure de meilleurs résultats de convergence théorique, ainsi qu'une meilleure performance en pratique.

La figure 2.3 illustre le fonctionnement simplifié d'un problème d'optimisation boîte noire. Des points d'essai $x \in \mathbb{R}^n$ sont fournis à la boîte noire par un logiciel d'optimisation. La boîte noire est évaluée et la valeur de la fonction objectif $f(x)$ est retournée au solveur, qui propose de nouvelles valeurs pour $x$.

<figure><img src="image_9.png" /><figcaption>Figure 2.3 Optimisation d’une boîte noire</figcaption></figure>

Souvent, les fonctions objectif des problèmes à optimiser sont très coûteuses à évaluer. Une fonction substitut simulant le même comportement que la vraie fonction est alors créée car elle est moins coûteuse à évaluer et devient la fonction à optimiser. Par exemple, cette façon de faire a été appliquée dans [4], afin de déterminer la localisation de capteurs de neige pour minimiser l'erreur d'interpolation de krigeage. Une bonne approximation des réserves de neige, surtout lors de la crue du printemps, est essentielle en production hydroélectrique afin d'éviter les déversements. La vraie fonction objectif de ce problème n'était pas disponible à l'utilisateur et c'est ce qui a motivé le choix d'une fonction substitut.

## Traitement des contraintes

Une façon de traiter les contraintes dans un problème d'optimisation boîte noire est la méthode de la barrière extrême. Pour un problème de minimisation, avec $f(x) : \mathbb{R}^n \to \mathbb{R}$ et
l'ensemble réalisable Ω, les contraintes sont traitées ainsi :

$$f(x) = \begin{cases} f(x) & \text{si } x \in \Omega \\ +\infty & \text{sinon.} \end{cases} \qquad (2.9)$$

Cette façon de traiter les contraintes implique de ne pas évaluer la boîte noire lorsqu'une des contraintes du problème est violée. D'autres méthodes existent, telles la méthode de la barrière progressive [11] ou du filtre [9], mais ne sont pas exposées dans ce document puisque la méthode de la barrière extrême est suffisante pour l'utilisation de l'optimisation boîte noire dans cette thèse.

### 2.3.4 Applications des méthodes d'optimisation boîte noire

Les méthodes d'optimisation boîte noire ont été appliquées à de nombreux problèmes d'ingénierie avec succès. Dans le domaine de l'hydrologie, l'optimisation boîte noire est utilisée par les auteurs de [4] afin de trouver la localisation optimale pour des GMON, un équipement utilisé pour mesurer l'équivalent en eau de la couverture de neige situé dans des endroits éloignés des bassins hydrographiques. Une autre étude [73] utilise l'optimisation boîte noire pour calibrer les 23 paramètres du modèle hydrologique HSAMI. Ce modèle est utilisé quotidiennement pour générer les prévisions d'apports pour les réservoirs. Plus spécifiquement appliqué aux turbines d'un système de production hydroélectrique, les auteurs de [5] utilisent l'optimisation boîte noire pour déterminer le nombre de turbines ainsi que les périodes exactes durant lesquelles elles peuvent être arrêtées afin d'être entretenues. Le troisième objectif de cette thèse consiste à optimiser les paramètres de la méthode de génération d'arbres de scénarios. Dans le même ordre d'idées, ces auteurs [12] utilisent l'optimisation boîte noire pour ajuster les valeurs des paramètres d'un algorithme.

De multiples applications existent dans d'autres domaines. Par exemple, les auteurs de [36] se servent de l'optimisation boîte noire pour déterminer des positions de semi-conducteurs sur un micro-processeur. Des fonctions de puissance différentes sont considérées et l'optimisation permet de minimiser la variation de température des semi-conducteurs selon les fonctions de puissance. Un problème non linéaire mixte en nombres entiers est résolu avec NOMAD [67]. Un design de moteur synchrone à aimants permanents doit minimiser le poids du moteur tout en maximisant le couple pour les vitesses de base et maximale. Plusieurs contraintes non linéaires sont considérées et les variables sont continues, comme le diamètre du stator ou entières, comme la grosseur du fil utilisé. D'autres applications en génie mécanique utilisent l'optimisation boîte noire pour optimiser la forme des pales d'une turbine [14] ou d'un hélicoptère [21]. Un compte-rendu [8] détaille l'optimisation boîte noire et de nombreux exemples
d'application sont présentés. Tous les exemples d'applications démontrent que l'optimisation
boîte noire permet de traiter des problèmes d'ingénierie complexes pour lesquels il serait
impossible, avec des méthode d'optimisation classiques telles l'optimisation linéaire ou l'op-
timisation en nombres entiers, de trouver une solution puisqu'ils n'ont pas de représentation
analytique des fonctions objectif ou des contraintes.
# CHAPITRE 3 ORGANISATION DE LA THÈSE

Le but de cette thèse est de développer un modèle d'optimisation stochastique à court terme pour le réseau hydroélectrique du Saguenay Lac-St-Jean. Dans le chapitre précédent, les méthodes d'optimisation existantes ont été présentées. Souvent, dans la littérature, des modèles déterministes sont utilisés, car les prévisions d'apports sont mises à jour régulièrement. Or, dans ce contexte, Rio Tinto souhaite utiliser un modèle stochastique, car des situations de prévisions qui ne se produisent finalement pas peuvent faire en sorte que les turbines soient forcées dans des conditions non efficaces, en plus de diminuer les niveaux de réservoirs inutilement. Dans un contexte opérationnel, le temps de calcul est primordial, alors la complexité des arbres de scénarios, qui servent à représenter l'incertitude dans le modèle d'optimisation à court terme, est étudiée.

En premier lieu, un modèle d'optimisation à court terme déterministe est développé. Les apports naturels des réservoirs sont supposés parfaitement connus. Pour ce faire, un algorithme de programmation dynamique est utilisé, avant toute optimisation, afin de créer des surfaces de puissance maximale pour chaque nombre de turbines en marche et chaque centrale. Ces surfaces sont ensuite approximées afin d'être utilisées comme fonction objectif du modèle d'optimisation. Un processus d'optimisation en deux phases est créé. La première phase résout le problème de chargement à l'aide d'un modèle non linéaire mixte en nombres entiers. La solution obtenue, soit le volume des réservoirs, le débit turbiné et le nombre de turbines en marche pour chaque centrale sont utilisés comme paramètres d'entrée de la deuxième phase. C'est un modèle linéaire en nombres entiers qui permet de déterminer la combinaison exacte de turbines en marche pour chaque centrale. Les résultats sont comparés aux décisions opérationnelles historiques, ainsi qu'aux solutions produites par un modèle mathématique en une phase, englobant le problème de chargement et de répartition. Ces travaux ont été publiés dans la revue *IEEE Transactions on Power Systems* [91] et se trouvent au chapitre 4.

En second lieu, le modèle d'optimisation déterministe est adapté afin de considérer l'incertitude liée aux apports naturels des réservoirs. Un modèle à court terme stochastique est développé. Une méthode de génération d'arbres de scénarios est retenue afin de représenter l'incertitude dans le modèle d'optimisation à court terme. Afin de valider les résultats, un banc d'essai est développé. Les prévisions d'apports sont mise à jour quotidiennement et un arbre de scénarios est généré. Le modèle d'optimisation stochastique est résolu, et seulement la solution du premier noeud dans l'arbre de scénarios est implémentée. Lorsque la vraie réalisation de l'apport est connue, les volumes des réservoirs sont mis à jour et le processus
est répété pour le nombre de jours désirés dans l'horizon de planification. Afin de valider
l'utilisation d'un modèle stochastique, les résultats sont comparés à l'utilisation du scéná-
rio médian des apports seulement dans le banc d'essai. Des essais avec plusieurs nombres
d'étapes dans l'arbre de scénarios ainsi que différents nombres de scénarios sont réalisés. Les
résultats sont présentés dans un article soumis à *European Journal of Operational Research*
[88] et se trouvent au chapitre 5.

En troisième lieu, la complexité des arbres de scénarios est étudiée. Le temps de calcul est un facteur important dans un contexte opérationnel comme celui-ci, et il est important d'obtenir une solution au modèle d'optimisation à court terme dans un délai raisonnable. Puisque seulement la solution du premier noeud de l'arbre est conservée, une comparaison entre l'utilisation d'un arbre de scénarios comme représentation de scénarios d'apports versus l'utilisation d'un peigne de scénarios est effectuée. Le même banc d'essai que pour le deuxième objectif est utilisé, et l'optimisation boîte noire est utilisée pour optimiser les paramètres de la méthode de génération d'arbres de scénarios afin de maximiser l'énergie produite. Une comparaison est aussi effectuée entre l'utilisation d'arbres de scénarios et l'utilisation du scénario médian dans le modèle d'optimisation à court terme. Les travaux sont publiés comme rapport technique dans la série des Cahiers du GERAD [90] ainsi qu'à un journal avec comité d'arbitrage et sont au chapitre 6.

Finalement, la synthèse des trois articles est présentée au chapitre 7 et la conclusion au
chapitre 8.
# CHAPITRE 4 ARTICLE 1 : SELF-SCHEDULING SHORT-TERM UNIT COMMITMENT AND LOADING PROBLEM

Cet article est publié :

S. Séguin, P. Côté et C. Audet, Self-scheduling short-term unit commitment and loading problem, *IEEE Transactions on Power systems*, 31 :1, 133-142, 2016.

**Abstract** - This paper presents a new method for solving the short-term unit commitment and loading problem for a specific hydropower system. Dynamic programming is used to compute maximum power output generated by a power plant. This information is then used as input of a two-phase optimization process. The first phase solves the relaxation of a nonlinear mixed-integer program in order to obtain the water discharge, reservoir volume and number of units working at each period in the planning horizon. The second stage solves a linear integer problem to determine which combination of turbines to use at each period. The goal is to maximize total energy produced over all periods of the planning horizon which consists of a week divided in hourly periods. Start-up of turbines are penalized. Numerical experiments are conducted on thirty different test cases for two Rio Tinto Alcan power plants with five turbines each.

**Index terms** - Hydro unit commitment and loading problem, optimization, non-linear programming, linear integer programming.

## 4.1 Introduction

The planning of hydroelectric systems is complex and requires different optimization processes. A good planning allows to produce more energy with the same quantity of water, generating substantial savings for the producer, even with a slight computational improvement [64]. Long-term, medium-term and short-term optimization models are used in order to manage the resources. Long-term optimization usually follows a few years of planning horizon [27],[71] and establishes future production potential under highly uncertain inflows in the basins. Medium-term optimization [47] is used to plan the reservoir volumes by estimating the quantity of water available for hydroelectric production on a weekly basis. Short-term optimization [40] is mandatory to determine how to split the available water volume in an optimal way between the turbines of a plant. Each turbine has a different efficiency curve
which means that for the same water discharge the power will differ. The planning horizon is a week divided in hourly periods and the problem consists of finding the water discharge as well as the volume of the reservoir for each plant in order to maximize energy production and penalize start-ups of turbines. The present work focuses on short-term optimization.

These optimization problems are difficult to solve since the hydroelectric production functions are nonconvex. They are also highly nonlinear and depend on turbine efficiency, net water head that is a nonlinear function of the water discharge and reservoir elevation and finally, water discharge of each unit. Furthermore, turbines have forbidden operating zones, which complicates the problem.

Short-term unit commitment and loading problem have been studied in the past and many researches are still undergoing. Many methods have been proposed to solve this problem, including dynamic programming. In [7], a power generation loss function is used to take into account the hydro power efficiency in the dispatch of the plants, but considers the forebay elevation insignificant and neglects the water balance equations. Another approach [109] maximizes basinwide operating efficiency and changes the scheduling of the units only when the energy generation or the water discharges exceed a deadband. These assumptions will cause the optimization solution to differ from the real ones. Formulations of interest must consider tailrace elevations, penstock losses as well as efficiency of the turbines, which are often set aside to simplify the problem. A different manner [23] approximates the influence of the head effect to linearize the power production function, while Ohishi [77] assumes that all units of a hydropower plant are similar and does not consider different power outputs being produced by the units. Another widely used method is the lagrangian relaxation to separate the linking constraints and solve subproblems that are easier to compute [42]. However, this method usually causes solutions to be slightly infeasible since the linking constraints are rarely satisfied with the first solutions. Heuristics are then used to obtain a feasible solution. Not only are the methods different, but the choice of the objective function can vary from one method to the other. Another formulation of the problem [59] is to minimize the sum of power losses and solve a relaxed mixed-integer nonlinear problem. Then, a simulation phase is processed to obtain a feasible solution for the relaxed constraints. Once again, a two-phase approach is necessary to obtain a feasible solution. In [97], the plant efficiency is maximized since it is known that water is not used in an efficient way to meet the demand in energy. Nonlinear approaches [28] have been considered, by linearizing the hydro power efficiency of the plants and the water level functions, but do not consider the unit commitment of the plants. Other techniques have been proposed, including genetic algorithms [79], ant algorithms [70] and network flows [78], but they require parameter tuning before obtaining solutions, which is not an issue when using mixed-integer formulations, nonlinear problems
or lagrangian relaxations, just to name a few. However these methods offer the advantage of treating efficiently nonconvex, nonsmooth and discontinuous problems.

This paper presents a new approach for modeling the short-term unit commitment and loading problem that requires a two-phase approach and allows to find a feasible solution at the end of the first stage. The first optimization dispatches generation among plants and seeks to maximize total energy production. The second stage uses this solution to select the unit commitment that maximizes total energy production but also penalizes unit start-ups. The models are then tested on two hydroelectric plants in series.

The paper is organized as follows. Section 4.1.1 presents the hydroelectric system studied in this paper and section 4.1.2 presents notation used throughout the paper. Section 4.2 gives an overview of the problem characteristics. Sections 4.3 presents the mathematical models developed which are a nonlinear model and a linear integer program. They respectively aim to solve the loading problem as well as the unit commitment problem. Extensive numerical experiments are reported in Section 4.4 on thirty test cases for two power plants with five turbines each, and concluding remarks are drawn in Section 4.5.

### 4.1.1 Saguenay-Lac-St-Jean hydroelectric system

The models presented in this paper are tested on the Saguenay-Lac-St-Jean hydroelectric system. It is privately owned by Rio Tinto Alcan in the province of Quebec. In this province, Hydro-Quebec is a government-owned corporation that generates, transmits and distributes electricity. Producers are reglemented and need to transit through Hydro-Quebec, which means that Rio Tinto Alcan cannot bid on the spot market. They sign contracts with the corporation that are limited in power and in energy and the same fares apply for their duration. This company operates aluminum plants in that region and can produce 90% of the energy they need to operate them. The installed capacity is of 3100 MW and is composed of 42 turbines divided in five hydroelectric plants situated on the Péribonka and Saguenay rivers. Five reservoirs are available and three of them have a stocking capacity of over 2000 $hm^3$. The hydrographic basins cover an area of about 75 000 $km^2$. For the remaining 10% of energy needs, energy is bought at a known fare from Hydro-Quebec. Historically, energy limits stated in the contracts are practically never reached and this is why the main purpose of the optimization is to produce the most energy possible with the available water. For the purpose of this paper, the models are tested on Chute-du-Diable and Chute-Savane plants that are both composed of five turbines. Specific constraints related to the system need to be taken into account when developing the model.
4.1.2 Notation

The following notation is used throughout the paper :

<table>
  <tr>
    <td>k &isin; {1, 2, ..., K}</td>
    <td>index of periods</td>
  </tr>
  <tr>
    <td>c &isin; {1, 2, ..., C}</td>
    <td>index of hydroelectric plants</td>
  </tr>
  <tr>
    <td>r &isin; {1, 2, ..., u<sup>c</sup>}</td>
    <td>index of hydroelectric plants upstream of plant c</td>
  </tr>
  <tr>
    <td>s &isin; {1, 2, ..., n<sub>k</sub><sup>c</sup>}</td>
    <td>index of surfaces corresponding to number of active turbines associated to hydroelectric plant c and period k</td>
  </tr>
  <tr>
    <td>l &isin; {1, 2, ..., n<sub>k</sub><sup>c</sup>}</td>
    <td>index of combinations associated to hydroelectric plant c and period k</td>
  </tr>
  <tr>
    <td>t &isin; {1, 2, ..., T<sup>c</sup>}</td>
    <td>index of turbines of hydroelectric plant c</td>
  </tr>
  <tr>
    <td>v<sub>k</sub><sup>c</sup></td>
    <td>volume of plant reservoir c at period k (hm<sup>3</sup>)</td>
  </tr>
  <tr>
    <td>q<sub>k</sub><sup>c</sup></td>
    <td>water discharge at plant c and period k (m<sup>3</sup>/s)</td>
  </tr>
  <tr>
    <td>g<sub>k</sub><sup>c</sup></td>
    <td>water spillage at plant c and period k (m<sup>3</sup>/s)</td>
  </tr>
  <tr>
    <td>&theta;</td>
    <td>start-up penalty for any turbine (MW)</td>
  </tr>
  <tr>
    <td>&beta;<sub>lk</sub><sup>c</sup></td>
    <td>power generated by combination l &isin; n<sub>k</sub><sup>c</sup> at plant c and period k</td>
  </tr>
  <tr>
    <td>y<sub>sk</sub><sup>c</sup> =</td>
    <td>
      <div style="display: flex; align-items: center;">
        <span style="font-size: 2.5em; line-height: 1;">{</span>
        <div style="margin-left: 1em;">
          <p>1 if surface s is chosen at period k</p>
          <p>for plant c</p>
          <p>0 otherwise</p>
        </div>
      </div>
    </td>
  </tr>
  <tr>
    <td>f<sub>lkt</sub><sup>c</sup> =</td>
    <td>
      <div style="display: flex; align-items: center;">
        <span style="font-size: 2.5em; line-height: 1;">{</span>
        <div style="margin-left: 1em;">
          <p>1 if turbine t of combination l</p>
          <p>for plant c is working at period k</p>
          <p>0 otherwise</p>
        </div>
      </div>
    </td>
  </tr>
  <tr>
    <td>x<sub>lk</sub><sup>c</sup> =</td>
    <td>
      <div style="display: flex; align-items: center;">
        <span style="font-size: 2.5em; line-height: 1;">{</span>
        <div style="margin-left: 1em;">
          <p>1 if combination l of plant c</p>
          <p>is chosen at period k</p>
          <p>0 otherwise</p>
        </div>
      </div>
    </td>
  </tr>
  <tr>
    <td>d<sub>tk</sub><sup>c</sup> =</td>
    <td>
      <div style="display: flex; align-items: center;">
        <span style="font-size: 2.5em; line-height: 1;">{</span>
        <div style="margin-left: 1em;">
          <p>1 if turbine t of plant c is started</p>
          <p>at period k</p>
          <p>0 otherwise</p>
        </div>
      </div>
    </td>
  </tr>
  <tr>
    <td>&chi;<sub>sk</sub><sup>c</sup>(v<sub>k</sub><sup>c</sup>, q<sub>k</sub><sup>c</sup>)</td>
    <td>power output function for surface s at period k and plant c (MW)</td>
  </tr>
</table>
<table>
  <tr>
    <td>δ<sup>c</sup><sub>k</sub></td>
    <td>inflow of plant c at period k (m<sup>3</sup>/s)</td>
  </tr>
  <tr>
    <td>w<sub>k</sub></td>
    <td>duration of period k (h)</td>
  </tr>
  <tr>
    <td>γ</td>
    <td>conversion factor from water discharge (m<sup>3</sup>/s) to (hm<sup>3</sup>/h)</td>
  </tr>
  <tr>
    <td>ζ<sub>k</sub></td>
    <td>conversion factor to energy units (GWh)</td>
  </tr>
  <tr>
    <td>v<sup>c</sup><sub>min</sub></td>
    <td>minimal volume of plant c reservoir (hm<sup>3</sup>)</td>
  </tr>
  <tr>
    <td>v<sup>c</sup><sub>max</sub></td>
    <td>maximum volume of plant c reservoir (hm<sup>3</sup>)</td>
  </tr>
  <tr>
    <td>q<sup>c</sup><sub>min</sub></td>
    <td>minimum water discharge at plant c (m<sup>3</sup>/s)</td>
  </tr>
  <tr>
    <td>q<sup>c</sup><sub>max</sub></td>
    <td>maximum water discharge at plant c (m<sup>3</sup>/s).</td>
  </tr>
</table>

## 4.2 Short-term unit commitment and loading problem

The short-term unit commitment and loading problems must determine a production plan to maximize energy production and penalize unit start-ups. This section describes the characteristics of the problem and presents an optimization model that requires a short computational time to solve.

### 4.2.1 Problem description

Power produced by a single hydroelectric generator [106] is given by the equation :

$$P = n(Q) \times g \times Q \times h_n \quad (4.1)$$

where *P* is the power output in *kW*, *n(Q)* is the turbine-generator overall efficiency, *g* is the gravitational acceleration in *m/s²*, *Q* is the turbine water discharge in *m³/s* and *h*<sub>n</sub> the net water head in *m*.

For a given turbine, power is a function of the water discharge, the net water head and the efficiency. Gross head is the difference between forebay and tailrace elevation. When water runs in the penstock, friction causes heat dissipation, thus diminishing energy. This phenomena causes a loss that needs to be considered in the power output calculation. Net water head is obtained from the gross head from which the losses are taken into account. Net water head is computed by :

$$h_n(v, Q) = h_f(v) - h_t(Q) - \psi(Q) \quad (4.2)$$

where *v* is the volume of the reservoir in *hm³*, *Q* is the total water discharge in *m³/s*, *h*<sub>f</sub> is a nonlinear function returning forebay elevation in *m*, *h*<sub>t</sub> is a nonlinear function returning
tailrace elevation in *m* and *ψ* is a nonlinear function returning friction losses in *m*.

The focus of this paper are Chute-du-Diable and Chute-Savane plants in the Saguenay-Lac-St-Jean hydroelectric system. A particularity of the system is that there is no analytical representation of the tailrace elevation for some plants. The value can only be calculated by a computer simulation. Since power depends on net water head and water discharge, there is no analytical representation of the power output functions. This particularity also deprives us of their derivatives.

Each turbine possesses its own efficiency curve, causing them to produce different power outputs for the same water discharge and net water head. Also, efficiency depends on water discharge of the turbine. Fig.4.1 illustrates the efficiency of a turbine as a function of water discharge for a given net water head at the Chute-du-Diable plant.

<figure><img src="image_10.png" /><figcaption>Figure 4.1 Turbine efficiency as a function of water discharge.</figcaption></figure>

In Fig.4.1, the water discharge of 48 *m³/s* generates the greatest power and is called maximum flow rate. This limit is fixed by the turbine manufacturer. When it is reached, water needs to be spilled, causing power and efficiency to decrease. Turbines should always be used within their efficiency zones for these reasons.

Another particularity that complicates the problem is the forbidden zones of operations of the turbines [50]. Under certain operating conditions, a vortex may occur in the turbine and create pressure variations that can damage components. These zones are forbidden and turbine are not operated when these conditions are met.
Also, unit restarts should be limited since they shorten equipment life [15]. There is also a cost associated to a unit start-up. It takes into account the history of expenses in maintenance and repairs in relation with the number of start-ups. This value is laboriously calculated by Rio Tinto Alcan and becomes a parameter in the optimization models.

One could formulate the energy production optimization problem as a linear integer program by discretizing unit water discharge, volumes and total water discharge for each turbine, power plant and period in order to maximize total energy production. However, the number of optimization variables would be extremely large. For instance, if total water discharge is discretized from 0 to 900 $m^3/s$, unit water discharge from 0 to 150 $m^3/s$ both with steps of 5 $m^3/s$, with 168 hourly periods for a week and volumes from 46 to 394 $hm^3$ discretized in 100 slices, then the number of binary variables would be of the order of $10^8$. Water discharges are discretized every $5 m^3/s$ since it is operationally impossible to obtain a finer precision, and volumes in 100 slices give a good final precision for the hydroelectric system. This suggests that the number of variables required is unrealistic for a real-time application.

### 4.2.2 Problem modeling

Power output of a single turbine is a function of two variables of the water discharge and the volume. However, there is a relation between the net water head and the volume of the reservoir. For the remainder of the paper, volume will be used to simplify notation. Total power output of a plant depends on total water discharge, number of working units and active turbines. Efficiency curves are specific for each turbine, hence unit power output is different for the same water discharge. Depending on the number of units working, but also on which units are employed, total power output is different. Instead of working directly with turbines in the model, fewer variables are needed if active turbines are grouped in combinations. For example, for the Chute-du-Diable plant, five turbines are available, but operational restrictions require a minimum of three active turbines. Table 4.2 lists the sixteen possible combinations. In each column, the numbers represent the actual active turbines.

<table><caption>Table 4.2 Turbine combinations at Chute-du-Diable</caption><thead><tr><th>3 active turbines</th><th>4 active turbines</th><th>5 active turbines</th></tr></thead><tbody><tr><td>123</td><td>145</td><td>1234</td></tr><tr><td>124</td><td>234</td><td>1235</td></tr><tr><td>125</td><td>235</td><td>1245</td></tr><tr><td>134</td><td>245</td><td>1345</td></tr><tr><td>135</td><td>345</td><td>2345</td></tr></tbody></table>

A model using combinations of active turbines needs to determine which one to use at each
period as well as the volume and the total water discharge. For a given volume and turbine, the power output depends on the water discharge. The same applies to combinations : if turbines in the combination differ, then the total power output will be different.

A dynamic programming algorithm is used to compute the combination power outputs. For every discretization of the water discharge, volume and combinations, the algorithm calculates power outputs. For a given power plant, volume and turbine, the states are $s^j = \{1, 2, ..., r\}$ remaining water to dispatch for turbine $j$ and the variables are $q^j = \{1, 2, ..., \min\{\bar{q}^j, Q\}\}$, unit water discharge for turbine $j$, where $Q$ is total water discharge and $\bar{q}^j$ is maximum unit water discharge. The objective is to maximize total power output. The solution is found recursively by solving a sub-problem for each turbine. For $j = n$, where $n$ is total number of turbines, optimal power output $f^{*j}(s^j) = P(s^j, v)$, where $v$ is the volume and $P$ is power output. The optimal unit water discharge is $q^{*j} = s^j$. For the other subproblems, that is for $j = n-1, n-2, ..., 1$, optimal power output $f^{*j}(s^j) = \max_{q^j} P(s^j, v) + f^{*j+1}(s^j - q^j)$ and optimal unit water discharge is $q^{*j} = s^j$ that maximizes $f^{*j}(s^j)$.

<figure><img src="image_11.png" /><figcaption>Figure 4.2 Power output at Chute-du-Diable.</figcaption></figure>

Fig.4.2a shows curves of the power output depending on the water discharge for all possible combinations at Chute-du-Diable power plant, for a volume of 376 $hm^3$ and Fig.4.2b is a close-up.

Observe that the power output decreases when the maximum efficiency of the turbine combination has been reached. These curves could be used to model the problem. The objective of the problem is to maximize total energy generation and a new function corresponding to the maximum of all combinations can be created. This new function returns the maximum power that can be generated for a given water discharge, a reservoir volume and a given number of units. It is represented by the bold curve on Fig.4.2a. This can be generalized for
the discretization of volumes. One hundred discretizations are done, between the minimum and maximum reservoir volume and are shown in Fig.4.3a. The interest for this function is that it gives us an upper bound on the optimal value of the problem for a volume and a given water discharge. The changes in the number of active turbines in the combination cause

<figure><img src="image_12.png" /><figcaption>Figure 4.3 Power output at Chute-du-Diable.</figcaption></figure>

a problem for modeling this surface since they correspond to nondifferentiable zones. The contour plot shown in Fig.4.3b illustrates this property. A way to overcome these difficulties is to create a surface for each number of active turbines in the combination. In this case, surfaces with three turbines working, four turbines and five turbines are created, as shown in Fig.4.4. Notice that in Fig.4.2a, the curves between 0% and 55% correspond to combinations with three active turbines, from 0% to 75% with four active turbines and five active turbines from 75% to 90%. As mentioned previously, a dynamic programming algorithm was used to compute power outputs for each water discharge, combination and volume discretization. Then, the maximum power output for every water discharge, volume and for each number of active turbines was retained to create the maximum power output surfaces for each number of active turbines. Smoothing splines were used to fit the power output data on the discretization points. The mathematical model must determine one surface per period, giving us at the same time the number of active turbines. Hence, if a turbine is unavailable for a given period, surfaces need to be recalculated without considering the unavailable turbine. Thus, for each plant and each period, the number of surfaces may vary, as well as the surfaces themselves. The dynamic programming algorithm is computed for every number of active turbines in the combination as well as for every combination of unavailable turbines and
<figure><img src="image_13.png" /><figcaption>Figure 4.4 Power output with different number of active turbines.</figcaption></figure>

power outputs are stored in a database. Models then consult the database, depending on
possible number of active turbines and available turbines. The volume of available water for
production is obtained from the medium-term optimization and is an input to the loading
problem, as well as the initial combination of working turbines. It is denoted $\hat{y}_{sk}^c$ and is a
parameter for the short-term unit commitment and loading optimization models. The use of
theses surfaces does not allow us to penalize the start-up of turbines, hence the optimization
needs to use a two-stage approach. The first stage is a nonlinear mixed-integer model and
returns the water discharge, reservoir volume as well as the number of active turbines. The
second stage is a linear integer model that determines exactly the combination of turbines to
use at each period in order to minimize the start-up of turbine so that energy is maximized.

A relevant contribution in this method is the pre-processing of the power output surfaces
before any optimization is conducted. This allows to decrease drastically the number of
optimization variables. Also, the creation of the power output surfaces results in smooth
surfaces that can be used in an optimization model. Therefore, derivatives are available and
these surfaces give us an analytical representation of the hydropower functions.

**4.3 Mathematical models**

This section presents the mathematical models used to solve the short-term unit commitment
and loading problem. The first one distributes generation among plants and the second
determines the active turbines.
### 4.3.1 Loading problem

The first optimization consists of a nonlinear mixed-integer program that determines the water discharge, reservoir volume and number of active turbines at each period for each power plant. The functions $\chi_{sk}^c(v_k^c, q_k^c)$ correspond to maximum power output of a given number of active turbines, for a given water discharge, volume, plant and period. These functions are smoothing splines that fit the data obtained from the dynamic programming algorithm. There is one function for each number of active turbines and period since turbines are not necessarily available all the time. The initial number of active turbines $\hat{y}_{s1}^c$ is known. Refer to Section 4.2.2 for more information on the calculation of maximum power output. The number of surfaces per plant and period is given by $n_k^c$.

Maximize total energy production :

$$ \max \sum_{c \in C} \sum_{k \in K} \sum_{s=1}^{n_k^c} \chi_{sk}^c(v_k^c, q_k^c) y_{sk}^c \zeta_k \quad (4.3) $$

subject to :

$$ \delta_k^c = v_{k+1}^c - v_k^c + \gamma w_k (q_k^c + g_k^c) \\ - \sum_{r=1}^{u^c} \gamma w_k (q_k^r + g_k^r) , \quad \forall c \in C, \forall k \in K \setminus \{1\} \quad (4.4) $$

$$ \sum_{s=1}^{n_k^c} y_{sk}^c = 1 , \quad \forall c \in C, \forall k \in K \setminus \{1\} \quad (4.5) $$

$$ y_{s1}^c = \hat{y}_{s1}^c , \quad \forall c \in C, \forall s = 1, 2, \dots, n_k^c \quad (4.6) $$

$$ v_{min}^c \le v_k^c \le v_{max}^c , \quad \forall c \in C, \forall k \in K \quad (4.7) $$

$$ q_{min}^c \le q_k^c \le q_{max}^c , \quad \forall c \in C, \forall k \in K \quad (4.8) $$

$$ q_k^c \ge 0 , \quad \forall c \in C, \forall k \in K \quad (4.9) $$

$$ v_k^c \ge 0 , \quad \forall c \in C, \forall k \in K \quad (4.10) $$

$$ g_k^c \ge 0 , \quad \forall c \in C, \forall k \in K \quad (4.11) $$

$$ y_{sk}^c \in B , \quad \forall c \in C, \forall k \in K, \\ \quad , \quad \forall s = 1, 2, \dots, n_k^c \quad (4.12) $$

$$ q_k^c, v_k^c, g_k^c \in R , \quad \forall c \in C, \forall k \in K. \quad (4.13) $$

Constraints (4.4) assure that water balance of the plants are met. The values of $v_1^c$ and $v_{k+1}^c$ are known and obtained from the medium-term optimization model. Constraints (4.5) force the model to choose only one surface at each period for each plant and constraints (4.6) feed the model with the turbines already working at the beginning of the planning horizon.
Constraints (4.7) and (4.8) represent the physical limits for volume of the reservoirs as well as the water discharge operational limits. Non-negativity of the variables are taken into account by constraints (4.9), (4.10) and (4.11). Finally, (4.12) impose binary variables and (4.13) real variables.

In practice, nonlinear mixed-integer programs require a large amount of computing time when solving since some variables, in this case those associated to the surfaces, are integer variables. Luckily, we can prove that solving the continuous relaxation of the surface variables $y_{sk}^c$ of this problem is sufficient to obtain an integer solution on these variables. Solving the relaxed nonlinear problem will return a solution with integer variables, even though imposing integer variables has been omitted.

**Proposition 4.3.1.** Solving the relaxation of integer program (4.3)-(4.13) leads to an integer solution for variables $y_{sk}^c$.

*Proof.* The proof of the result is done by showing that the matrix of constraints for variables $y_{sk}^c$ is totally unimodular. Problem (4.3)-(4.13) can be re-written as follows :

$$ \max \sum_{\omega \in \Omega} \sum_{c \in C} \sum_{k \in K} \sum_{s=1}^{n_k^c} \eta_{\omega k} y_{sk}^c \zeta_k \quad (4.14) $$

subject to :

$$ \sum_{s=1}^{n_k^c} y_{sk}^c = 1 \quad , \quad \forall c \in C, \forall k \in K \setminus \{1\} \quad (4.15) $$

$$ y_{s1}^c = \hat{y}_{s1}^c \quad , \quad \forall c \in C, \forall s = 1, 2, \dots, n_k^c \quad (4.16) $$

$$ y_{sk}^c \in B \quad , \quad \forall c \in C, \forall k \in K, \forall s = 1, 2, \dots, n_k^c. \quad (4.17) $$

where $\eta_{wk}$ is the power for $\omega$ the feasible set with respect to constraints (4.4) and $\hat{y}_{sk}^c$ is the initial combination of turbines working from the medium-term optimization. Since the variables $y_{sk}^c$ do not appear in the water balance constraints, the model is written with the feasible set satisfying constraints (4.4). Thus, $\eta_{wk}$ is equivalent to $\chi_{sk}^c(v_k^c, q_k^c)$, but for simplification, $v_k^c$ and $q_k^c$ do not appear in the model to make it clear that variables $y_k^c$ only influence the choice of the surface constraints.

Denote the matrix of coefficients of the constraints for problem (4.14)-(4.17) by A. Wolsey [107] shows A is totally unimodular if and only if the following three conditions are satisfied :

1. $a_{ij} \in \{+1, -1, 0\} \forall i, j$.

2. Each column of A contains at most two nonzero coefficients ($\sum_{i=1}^m |a_{ij}| \le 2$).
3. There exists a partition $(M_1, M_2)$ of the set $M$ of rows of $A$ such that each column $j$ containing two nonzero coefficients satisfies $\sum_{i \in M_1} a_{ij} - \sum_{i \in M_2} a_{ij} = 0$.

Condition 1. is satisfied by Equations (4.15-4.16). Each column of the matrix $A$ has a single element, which imply that Conditions 2. and 3. are satisfied.

Therefore, $A$ is totally unimodular and there is an integer optimal solution of the continuous relaxation. $\square$

The following nonlinear relaxed program can be solved, where $Y_{sk}^c$ are continuous variable associated to the surfaces. Maximize total energy produced at each period :

$$ \max \sum_{c \in C} \sum_{k \in K} \sum_{s=1}^{n_k^c} \chi_{sk}^c (v_k^c, q_k^c) Y_{sk}^c \zeta_k \quad (4.18) $$

subject to :

$$ (4.4) \qquad (4.19) $$

$$ \sum_{s=1}^{n_k^c} Y_{sk}^c \le 1, \quad \forall c \in C, \forall k \in K \qquad (4.20) $$

$$ (4.6) - (4.10) \qquad (4.21) $$

$$ Y_{sk}^c \ge 0, \quad \forall c \in C, \forall k \in K, \forall s = 1, 2, \dots, n_k^c \qquad (4.22) $$

$$ Y_{sk}^c \in \mathbb{R}, \quad \forall c \in C, \forall k \in K, \forall s = 1, 2, \dots, n_k^c. \qquad (4.23) $$

Constraints remain the same, except for (4.20) that becomes an inequality and (4.22) and (4.23) that are the non-negativity constraints for continuous variables associated to the choice of the surface.

### 4.3.2 Unit commitment

The solution produced by the nonlinear relaxed program is an input to the unit commitment problem. The unit commitment model is a linear integer program and it determines the exact combination of turbines to use in order to maximize total energy production at each period and penalize start-up of turbines. The initial combination of working turbines $\hat{x}_{lk}^c$ is known. From the loading problem, water discharge, volume and number of active turbines are obtained for each period and plant. Power output is then computed with a dynamic programming algorithm for every combination corresponding to the number of active turbines, for given values of volume and water discharge calculated from the loading problem. These power output values are stored in $\beta_{lk}^c$, which are parameters in the unit commitment model.
The number of combinations for a given period and power plant is given by $n_k^c$ and number of turbines of a hydroelectric plan by $T^c$.

The optimization problem maximizes energy produced and penalizes turbine start-ups :

$$ \max \sum_{c \in C} \sum_{k \in K} \sum_{l=1}^{n_k^c} \beta_{lk}^c x_{lk}^c \zeta_k - \sum_{c \in C} \sum_{k \in K} \sum_{t=1}^{T^c} d_{tk}^c \theta \zeta_k \quad (4.24) $$

subject to :

$$ \sum_{l=1}^{n_k^c} x_{lk}^c = 1, \quad \forall c \in C, \forall k \in K \setminus \{1\} \quad (4.25) $$

$$ x_{lk}^c f_{lkt}^c - x_{lk-1}^c f_{lk-1t}^c \le d_{tk}^c, \quad \forall c \in C, \forall k \in K, \forall l = 1, 2, \dots, n_k^c \\ , \quad \forall t = 1, 2, \dots, T^c \quad (4.26) $$

$$ x_{l1}^c = \hat{x}_{l1}^c, \quad \forall c \in C, \forall l = 1, 2, \dots, n_k^c \quad (4.27) $$

$$ x_{lk}^c \in B, \quad \forall c \in C, \forall k \in K, \\ , \quad \forall l = 1, 2, \dots, n_k^c \quad (4.28) $$

$$ d_{kt}^c \in B, \quad \forall c \in C, \forall k \in K, \\ , \quad \forall l = 1, 2, \dots, n_k^c \\ , \quad \forall t = 1, 2, \dots, T^c. \quad (4.29) $$

The constraints (4.25) ensure that only one combination is chosen at each period. Constraints (4.26) are the linking constraints between start-up variables and combination choice. Constraints (4.27) force the initial combination of turbines working. Finally, constraints (4.28) and (4.29) are the declaration of binary variables.

One of the advantages of solving separately the unit commitment problem is that there is no approximation of the data, since $\beta_{lk}^c$ is calculated with a dynamic programming algorithm. Also, the advantage of a two-phase optimization approach is that it is much easier to solve a continuous nonlinear model followed with a linear integer model than a single nonlinear mixed integer problem.

## 4.4 Computational results

The mathematical models of Sections 4.3.1 and 4.3.2 are tested on two of the five hydroelectric plants that compose the Saguenay-Lac-St-Jean hydroelectric system. These two plants, Chute-du-Diable and Chute-Savane, are in series and both have five turbines. They were chosen since they are the smallest sub-system with two plants in series and the motivation is to validate the optimization developed before expanding to the whole hydroelectric system.
Also, the models developed are deterministic and the aim of subsequent studies will be to consider uncertainty in the weather forecasts.

The current procedure used at Rio Tinto Alcan to schedule hydro power production is a gradient-based method with linear constraints. A short-term optimization software computes total water discharge at each plant and period. For the unit commitment, a dynamic programming algorithm combined with the engineer's experience are used to make the final decision. The purpose of this paper is to develop a fully automatic procedure that gives better, or at least as good results as currently. The planning horizon of the models is partitioned into 168 hourly periods for one week. Thirty weekly scenarios are tested, all provided from the historical database. Since data is available every two minutes, a mean per hour of the water discharge and the volume of the reservoirs is calculated in order to compare results. Initial volume, final volume and initial combination of working turbines are provided to the model to make the best possible comparison. The same starting point, in the middle of volume and water discharge discretizations is given to the model as an initial solution. Inflows of the basins are assumed known and the computed model is deterministic.

The Ipopt [105] nonlinear solver is used for the first nonlinear program and Xpress-MP [2] for the second linear program. Numerical experiments are conducted with an Intel Xeon Processor E5-2650, with 8 cores at 2 GHz and 61 Go of RAM memory. The nonlinear program has a total of 1680 real variables and 674 constraints. For the second linear integer program, 7056 binary variables and 26918 constraints are necessary. The number of variables and constraints for the second model will be less if some turbines are unavailable due to maintenance or repair at certain periods. Numerical comparisons of our approach versus the real values show an improvement on 27 out of the 30 test cases, ranging from 0.002 GWh to 2.145 GWh. The average improvement over all cases is of 0.4 GWh. In the province of Quebec, a 1 MWh earning represents roughly 40$ in savings for the producer. In this particular case, an improvement of 0.4 GWh translates into 832 000$ savings for a year.

Our approach is slightly sensitive to the starting point value. We observed that for the three cases in which our approach did not improve the solution, a better solution was found by changing the starting point. The computational time to solve the unit commitment model is very low. The loading problem takes an average of 1.41 seconds and the longest time is 7.71 seconds. Results obtained with the same initial solution are shown in Table 4.3. For each hydroelectric plant, total energy, with penalties due to start-ups, both for solution obtained with the optimizer and the historical database are listed. Also, the difference between optimized solutions and the historical database are computed. A positive value indicates that the optimizer produces a better solution than reality and a negative value indicates the opposite. Also, the total difference of start-ups of turbines is listed. A positive difference
indicates optimized solution has more starts than real case and a negative difference indicates
real case has more starts.

<table><caption>Table 4.3 Total energy production with same initial solution</caption><thead><tr><th rowspan="2">#</th><th colspan="2">Chute-du-Diable</th><th colspan="2">Chute-Savane</th><th colspan="2">Total</th><th colspan="2">Diff.</th></tr><tr><th>Optimized<br>GWh</th><th>Real<br>GWh</th><th>Optimized<br>GWh</th><th>Real<br>GWh</th><th>Optimized<br>GWh</th><th>Real<br>GWh</th><th>GWh</th><th>#<br>start</th></tr></thead><tbody><tr><td>1</td><td>27,12</td><td>26,82</td><td>38,09</td><td>37,93</td><td>65,21</td><td>64,75</td><td>0,46</td><td>0</td></tr><tr><td>2</td><td>31,96</td><td>31,97</td><td>35,41</td><td>35,25</td><td>67,37</td><td>67,21</td><td>0,16</td><td>1</td></tr><tr><td>3</td><td>33,42</td><td>33,21</td><td>36,55</td><td>36,21</td><td>69,97</td><td>69,42</td><td>0,55</td><td>-1</td></tr><tr><td>4</td><td>22,52</td><td>22,07</td><td>24,73</td><td>23,80</td><td>47,25</td><td>45,87</td><td>1,38</td><td>0</td></tr><tr><td>5</td><td>25,01</td><td>25,05</td><td>28,31</td><td>28,12</td><td>53,32</td><td>53,17</td><td>0,15</td><td>-3</td></tr><tr><td>6</td><td>29,36</td><td>29,07</td><td>31,68</td><td>30,98</td><td>61,03</td><td>60,05</td><td>0,98</td><td>2</td></tr><tr><td>7</td><td>22,74</td><td>22,87</td><td>30,00</td><td>29,96</td><td>52,74</td><td>52,82</td><td>-0,08</td><td>5</td></tr><tr><td>8</td><td>32,26</td><td>31,56</td><td>30,18</td><td>29,67</td><td>62,44</td><td>61,23</td><td>1,21</td><td>3</td></tr><tr><td>9</td><td>27,29</td><td>26,78</td><td>38,38</td><td>38,42</td><td>65,67</td><td>65,20</td><td>0,47</td><td>1</td></tr><tr><td>10</td><td>32,06</td><td>31,91</td><td>36,10</td><td>35,52</td><td>68,16</td><td>67,42</td><td>0,74</td><td>0</td></tr><tr><td>11</td><td>36,67</td><td>36,68</td><td>38,41</td><td>38,39</td><td>75,07</td><td>75,07</td><td>0,00</td><td>0</td></tr><tr><td>12</td><td>28,68</td><td>28,30</td><td>36,59</td><td>36,14</td><td>65,27</td><td>64,44</td><td>0,83</td><td>0</td></tr><tr><td>13</td><td>33,94</td><td>33,93</td><td>38,80</td><td>38,72</td><td>72,73</td><td>72,65</td><td>0,08</td><td>0</td></tr><tr><td>14</td><td>26,52</td><td>26,15</td><td>36,76</td><td>36,59</td><td>63,28</td><td>62,74</td><td>0,54</td><td>1</td></tr><tr><td>15</td><td>24,04</td><td>23,54</td><td>35,29</td><td>35,35</td><td>59,33</td><td>58,88</td><td>0,45</td><td>1</td></tr><tr><td>16</td><td>33,80</td><td>33,53</td><td>37,85</td><td>37,63</td><td>71,66</td><td>71,16</td><td>0,50</td><td>0</td></tr><tr><td>17</td><td>31,13</td><td>31,13</td><td>36,54</td><td>36,43</td><td>67,67</td><td>67,56</td><td>0,11</td><td>0</td></tr><tr><td>18</td><td>30,18</td><td>29,89</td><td>31,84</td><td>31,94</td><td>62,01</td><td>61,83</td><td>0,18</td><td>2</td></tr><tr><td>19</td><td>29,10</td><td>28,70</td><td>31,26</td><td>31,37</td><td>60,36</td><td>60,06</td><td>0,29</td><td>4</td></tr><tr><td>20</td><td>28,57</td><td>28,26</td><td>31,33</td><td>31,32</td><td>59,90</td><td>59,58</td><td>0,32</td><td>3</td></tr><tr><td>21</td><td>27,04</td><td>26,96</td><td>29,80</td><td>29,71</td><td>56,84</td><td>56,67</td><td>0,17</td><td>1</td></tr><tr><td>22</td><td>28,17</td><td>27,83</td><td>30,68</td><td>30,27</td><td>58,84</td><td>58,10</td><td>0,75</td><td>-4</td></tr><tr><td>23</td><td>25,05</td><td>25,11</td><td>27,77</td><td>27,57</td><td>52,82</td><td>52,68</td><td>0,14</td><td>0</td></tr><tr><td>24</td><td>33,64</td><td>32,67</td><td>36,97</td><td>35,80</td><td>70,61</td><td>68,47</td><td>2,15</td><td>-1</td></tr><tr><td>25</td><td>20,54</td><td>20,64</td><td>22,90</td><td>22,78</td><td>43,44</td><td>43,42</td><td>0,02</td><td>2</td></tr><tr><td>26</td><td>31,91</td><td>31,87</td><td>35,14</td><td>35,14</td><td>67,06</td><td>67,01</td><td>0,05</td><td>-2</td></tr><tr><td>27</td><td>20,34</td><td>20,48</td><td>23,41</td><td>23,25</td><td>43,74</td><td>43,73</td><td>0,02</td><td>3</td></tr><tr><td>28</td><td>29,48</td><td>29,62</td><td>32,83</td><td>32,79</td><td>62,31</td><td>62,41</td><td>-0,11</td><td>-2</td></tr><tr><td>29</td><td>33,34</td><td>33,34</td><td>36,68</td><td>36,56</td><td>70,02</td><td>69,90</td><td>0,12</td><td>1</td></tr><tr><td>30</td><td>29,08</td><td>29,27</td><td>31,83</td><td>31,85</td><td>60,91</td><td>61,12</td><td>-0,21</td><td>3</td></tr></tbody></table>

Fig.4.5 is a histogram comparing energy gain between optimized solution and real test cases. In twelve out of the thirty test cases, the optimized solution improved the quantity of energy produced between 0 and 0.2 GWh. For ten other cases, the improvement is between 0.2 and 1.0 GWh. Three cases exceed 1.0 GWh and the largest gain is close to 2.2 GWh. In 27 of these cases the improvement is positive, and in the other 3 the improvement is slightly
negative.

<figure><img src="image_14.png" /><figcaption>Figure 4.5 Histogram of differences of energy.</figcaption></figure>

### 4.4.1 Interpretation of the results

We illustrate the differences between the real and the optimized solutions by analyzing two of the test cases. These cases were selected since one of them proposes a very different production plan than the decisions actually taken at the moment. Case 1 fills the reservoir during the week and case 6 has similar volumes at the beginning and the end of the week.

Since power is a two variable function of the net water head and water discharge, both graphs are presented. For Chute-du-Diable and Chute-Savane, a graph of the water discharge comparing optimized solution and reality as well as a graph comparing water heads are displayed for both cases.

Let us analyze case 1 in details. Fig.4.6a shows the results for Chute-du-Diable power plant. The graphic on top shows the total water discharge of the plant in $m^3/s$ and the graphic below shows the net water head in $m$. The period for both graphs is 2 minutes, as in the historic database. The optimized result over the entire week in hourly periods are transposed every 2 minutes to allow a visual comparison between historic and optimized results. Even though a mean per hour of the 2 minute inflows has been computed for the optimization, the total volume is equivalent for 2 minutes real results and 1 hour optimized solutions. For each graph, the optimized solution is presented with a dashed line and the real cases with a filled line. Fig.4.6b shows the results for Chute-Savane. The top graphic of Fig.4.6a shows the optimized solution in which the total water discharge at the plant is 670 $m^3/s$ during
the whole planning horizon. Graphic on the bottom, on Fig.4.6a, illustrates the water head is similar throughout all the week for database and optimizer solution. Fig.4.6b on top shows database solution give a constant water discharge of about 780 $m^3/s$ with a higher discharge at the beginning of the week, and that optimized solution has also a constant solution with a peak of the water discharge around period 1000. As for the graphic on the bottom, the optimized water head is higher than database results throughout the week. The solution provided by the optimizer produces an improvement of 0.462 GWh.

<figure><img src="image_15.png" /><figcaption>Figure 4.6 Case 1 water discharge and water head.</figcaption></figure>

For case 6, the solutions for both power plants have strategies different than decisions made at that moment. Chute-du-Diable's reservoir is lowered then filled and Chute-Savane is filled, then lowered and filled once again. Results are presented on Fig.4.7.

These two cases were selected to illustrate important differences. Case 1 fills the reservoir throughout the week and case 6 keeps reservoir at the same level. These graphs show that the optimized strategy differ from what was done in reality, and improve the production. Case 1 produces 0.462 GWh more and case 6, 0.982 GWh more. Our models allow us to produce a solution within a very satisfying computational time of a few seconds. Also, it demonstrates that the proposed optimization lead to a fully automatic procedure that allows to obtain results at least as good as the current procedure.

We solved the short-term unit commitment and loading problem with a deterministic model, which means that no uncertainties are taken into account on the inflows of the basin. It gives more liberty to the optimization to vary the reservoir volumes knowing exactly water inflows
<figure><img src="image_16.png" /><figcaption>Figure 4.7 Case 6 water discharge and water head.</figcaption></figure>

that will occur at the next period. Computational results show that the model proposed could be modified and extended to take into account the uncertainties related to the weather forecasts.

### 4.4.2 Benchmarking of the proposed solution

The current procedure at Rio Tinto Alcan combines a gradient-based method with the engineer's experience. Thus, to validate the proposed optimization method, it is necessary to compare with a mathematical algorithm that does not apply a post-processing of the results. In a perfect situation, the output of the gradient-based method could have been used, but the model does not perform well. Hence, the results obtained with our two-phase optimization model would not have been a fair comparison since it always performs much better.

We propose the following benchmark: the two-phase models are combined into a single a nonlinear mixed integer model and solved with the Bonmin [20] software package. The same dynamic programming algorithm as in Section 4.3.1 is used but instead of computing surfaces for every number of active turbines in the combination, they are computed for every combination of turbines since the optimization must be carried out in a single phase. The continuous variables are: water discharge, water spilled and volume at each period for each plant. The binary variables are: combination used and startup of turbines at each period, for each plant. The constraints are the same as for the two-phase models combined. The objective function is to maximize total energy production and to penalize startup of turbines.
The benchmark test cases are the ones from Section 4.4 in which all ten turbines are available at each period of the planning horizon. There are ten such test cases. Table 4.4 summarizes the results.

<table><caption>Table 4.4 Benchmarking on instances with 10 turbines</caption><thead><tr><th rowspan="2">Case</th><th colspan="2">Benchmark</th><th colspan="2">two-phase</th></tr><tr><th>GWh</th><th>Time m :s</th><th>GWh</th><th>Time m :s</th></tr></thead><tbody><tr><td>1</td><td>-</td><td>&gt; 12m</td><td>65.21</td><td>0 : 01</td></tr><tr><td>7</td><td>52.67</td><td>0 : 07</td><td>52.74</td><td>0 : 02</td></tr><tr><td>9</td><td>-</td><td>&gt; 12m</td><td>65.67</td><td>0 : 01</td></tr><tr><td>11</td><td>75.29</td><td>0 : 09</td><td>75.07</td><td>0 : 01</td></tr><tr><td>13</td><td>72.73</td><td>0 : 11</td><td>72.73</td><td>0 : 01</td></tr><tr><td>14</td><td>63.09</td><td>2 : 53</td><td>63.28</td><td>0 : 02</td></tr><tr><td>15</td><td>59.05</td><td>2 : 33</td><td>59.33</td><td>&lt; 1s</td></tr><tr><td>17</td><td>67.38</td><td>0 : 09</td><td>67.67</td><td>&lt; 1s</td></tr><tr><td>22</td><td>58.58</td><td>0 : 09</td><td>58.84</td><td>0 : 01</td></tr><tr><td>28</td><td>61.88</td><td>0 : 25</td><td>62.31</td><td>0 : 06</td></tr></tbody></table>

Out of the ten cases, two cannot be solved by Bonmin with the default parameters. After 12 minutes, the solver aborted with an infeasible solution, while the two-phase optimization process solves these both cases in one second. For two other cases, the benchmark required more than 2 minutes and again, the proposed two-phase model solves these two cases in less than 2 seconds. Finally, the proposed method in this paper solves the remaining cases in less than 6 seconds, which is inferior to the fastest time for Bonmin.

For the test cases solved by both the benchmark and two-phase methods, the final objective function values are comparable, as seen in Table 4.4.

One may notice the higher values of the objective function for the two-phase optimization process than for the benchmark. As said previously, the surfaces were computed for every combination of turbines in the benchmark model. Since all the optimization, which are the loading and the unit commitment model, needs to be carried out in a single phase, polynomial equations were used to fit the data. In the two-phase process, the first phase power surfaces were fitted with splines, but in the second phase, the value of power for a given water discharge and volume were computed using a dynamic programming algorithm, which means that there was no approximation on the data. Therefore, there is no way of overcoming this difficulty when solving the models in a single phase, so that explains the differences in the objective value. The polynomial equations used to fit the data in the benchmark will inevitably have over or under estimations, which is not the case with the two-phase optimization. Fig.4.8 illustrates the water heads for both plants and both algorithmic approaches. The results are
very similar.

<figure><img src="image_17.png" /><figcaption>Figure 4.8 Case 9 benchmark and two-phase approach water head.</figcaption></figure>

The comparison between the output of a nonlinear mixed integer model and the two-phase optimization method developed in this paper shows that our method solves all the proposed test cases in a computational time that is much faster than in a one phase optimization process. The computational times increase with the number of turbines. In our test cases, problems have 10 turbines, but we must keep in mind that some hydroelectric plants have up to 18 turbines, and the single phase approach with the Bonmin solver cannot solve all problems.

## 4.5 Conclusion

Short-term unit commitment and loading problems for the Saguenay-Lac-St-Jean system are complex to solve since the hydroelectric production functions are nonconvex and nonlinear and we do not have analytical representations of them. We have proposed a model with a reasonable number of variables, embedded into a two-stage optimization approach. The first stage solves the relaxation of a nonlinear mixed-integer program in order to find volume, water discharge and number of active turbines at each period. The second stage solves a linear integer model to find the exact combination of turbines that maximizes total energy but also penalizes start-up of turbines. Dynamic programming is used to calculate total power output that can be generated by a certain combination of active turbines. This data is then used as parameters for both models. The pre-processing of the power output surfaces is a
main contribution to this paper since it decreases the number of optimization variables. Also,
even though tailrace elevations have no analytical representation, the power output surfaces
allows us to use solvers that require derivatives, since they are smooth surfaces and they give
a good representation of the hydropower production functions. The approach proposed in this
paper allows us to find a solution in a computational time that is more than satisfying for
needs of operation and was compared to a large mixed integer nonlinear problem formulation.
Multistarts or variable neighborhood searches [53] will be the subject of future research. Also,
other developments based on this method will involve using uncertainty related to inflows in
order to create a stochastic programming model.

Acknowledgment

The authors would like to thank Marco Latraverse from Rio Tinto Alcan for his implication
and advice throughout this project. Also, the comments of the referees were very helpful and
improved the quality of the paper.
# CHAPITRE 5 ARTICLE 2 : STOCHASTIC SHORT-TERM HYDROPOWER PLANNING WITH INFLOW SCENARIO TREES

Cet article a été soumis à la revue *European Journal of Operational Research* le 29 septembre 2015. Une révision de l'article a été soumise le 19 avril 2016 et est toujours en processus de révision.

Le rapport technique suivant est publié :

S. Séguin, S.-E. Fleten, A. Pichler, P. Côté et C. Audet, Stochastic short-term hydropower planning with inflow scenario trees, *Cahiers du GERAD*, (G-2015-97), 2016, Révisé avril 2016.

**Abstract :** This paper presents an optimization approach to solve the short-term hydro-power unit commitment and loading problem with uncertain inflows. A scenario tree is built based on a forecasted fan of inflows, which is developed using the weather forecast and the historical weather realizations. The tree-building approach seeks to minimize the nested distance between the stochastic process of historical inflow data and the multistage stochastic process represented in the scenario tree. A two-phase multistage stochastic model is used to solve the problem. The proposed approach is tested on a 31 day rolling-horizon with daily forecasted inflows for three power plants situated in the province of Quebec, Canada, that belong to the company Rio Tinto.

**Key Words :** Large scale optimization, nonlinear programming, OR in energy, scenarios, stochastic programming.

**Acknowledgments :** The authors would like to thank Marco Latraverse, and Rio Tinto, for providing data necessary to this study. This work was supported by NSERC, FRQNT and Rio Tinto. Also, a grant for international mobility, awarded by FRQNT through GERAD allowed the research to be conducted at the Norwegian University of Science and Technology. Stein-Erik Fleten acknowledges financial support from the Research Council of Norway through project 243964/E20. Alois Pichler gratefully acknowledges support of the Research Council of Norway through grant 207690/E20. The authors are thankful to two anonymous referees, for providing constructive comments that greatly improved the paper.
## 5.1 Introduction

Hydroelectric producers invest time and resources in developing optimization tools to gain efficiency in the use of water, since even small improvements lead to significant savings. Short-term optimization is used at the power plant level to dispatch available water for production between the turbines. Each turbine has a different efficiency. The amount of water available for production, or reservoir trajectories, is determined from the medium-term optimization and considers demand, uncertainty in the inflows, and travel time of the water between the plants. Short-term optimization is often considered to be deterministic [99] by making the assumption that the inflows are known [41] or by neglecting water balance constraints [7] at such a short time scale, but does not allow planning under different forecasts. Also, [86] have shown that considering uncertainty in short-term decision models may lead to improvements.

The focus of this paper is stochastic optimization applied to the short-term hydropower optimization problem. By considering uncertain inflows, turbines will be used in a more efficient manner since the stochastic model results in a compromise between high and low forecasted inflows. For example, in situations where reservoirs are nearly full, considering uncertain inflows when high inflows are expected prevents lowering the reservoir and force turbines into inefficient zones, which results in energy production loss in the future if these high inflows do not occur.

Few papers have looked specifically into short-term hydropower models with uncertain inflows. In [91], a short-term hydropower optimization model treats deterministic inflows. Water head variations are considered and nonlinearities and nonconvexities of the hydropower production function are accounted for. In [45], uncertainty of prices and inflows is considered. The authors use time series analysis to model the water inflows, which is represented by a scenario tree in the stochastic programming model. Start-up costs are considered and a multistage stochastic model is approximated by a two-stage model. A mixed-integer linear program is used. The net water head is assumed to vary with the water discharge only, so hydropower production functions depend only on the water discharge.

In [83], the only uncertainty considered is demand. The deterministic model is a linear integer model, which is an approximation of a nonlinear mixed integer model. Once again, the hydro-power production function depends only on water discharge. For some hydropower systems, neglecting the water head is not a possible avenue since many of the reservoirs have small capacities. Consequently, the water head effect is important in a short-term optimization, even with short time steps.

Many assumptions are made when solving the short-term unit commitment model, since they
are complex to solve. They have a large amount of variables, power production functions are nonlinear and efficiency is different for every turbine. The most common assumption is to neglect water head variations leading to linear power production functions.

When uncertainty arises and one wants to solve the optimization models, two main streams of ideas have been applied in the optimization community. Stochastic dynamic programming has been used extensively to solve hydropower optimization models [95, 100], as well as variants such as sampling stochastic dynamic programming [32] or stochastic dual dynamic programming [92]. These models are well suited for long or medium-term horizons but for short-term models, the state space is huge and it is very difficult, if not impossible, to solve them. In order to prevent the optimization process to empty out the reservoirs in the short-term model, values are assigned to the remaining water at the end of the planning horizon, which are obtained with stochastic dynamic programming or stochastic dual dynamic programming for example. In [69], a new method to generate inflows, based on periodic autoregressive models, is used as input to a stochastic dual dynamic programming algorithm that allows to schedule a hydro-thermal system located in Brazil.

The other stream is stochastic programming. A two-stage stochastic model [19] consists of two stages of decisions. The first-stage decisions need to be taken without knowing the realization of the uncertainty in the future, while the second stage decisions are taken when the uncertainty is revealed.

Usually, uncertainty is represented by scenarios. Each scenario is a possible realization of the uncertainty. Multiple scenario generation methods have been used in the past to approximate the distributions of the stochastic parameters. An overview of these methods, as well as evaluating the quality of a scenario tree is found in [61]. In [34], a periodic autoregressive process is used to fit historical data of the prices and to generate prices for the stochastic model. The scenario tree is built by sampling the distribution fitted with the model for the different nodes. Another method creates a discrete distribution of the uncertain parameter by matching some specific statistical properties. In [56], the first four moments, mean, variance, skewness and kurtosis are matched. Multiple pitfalls arise from this method and one must ensure the scenario tree represents possible outcomes of the uncertainty. A survey of techniques for generating scenario trees appears in [35] and includes recombining of data paths, contamination method and matching. Also, copulas have been used to generate scenarios for two-stage stochastic problems [61]. This method offers the advantage of treating dependencies better than with correlation alone. Other methods are scenario reduction [38, 52]. An initial scenario tree is required and forward selection, or backward reduction is applied in order to reduce it and minimize the computational time required to solve the stochastic optimization
model. The effect of the reduction on the solution accuracy, applied to a cascaded system of hydropower reservoirs is found in [108].

Other methods to deal with uncertainty on the inflows include robust optimization techniques [13] and probabilistic constrained programming [103]. Robust optimization solves models that have uncertain parameters over uncertainty sets. Therefore, the optimization seeks to find a solution that is feasible regardless of the outcome of the uncertainty. In [6], a rolling-horizon scheme is used and robust optimization is applied to the decision of day 1 while the rest of the horizon is considered deterministic. This is interesting as the uncertainty is applied to the important decisions. A drawback of robust optimization is the formulation of the uncertainty. In the historical records, some values of inflows may be very low and others very high. Therefore, it is difficult to define what are the best bounds for the uncertainty set, as well as capturing any nonlinear dynamics present. In probabilistic constrained programming, constraints are to be respected given a certain probability. A cascaded hydropower system is optimized with probabilistic constrained programming in [103]. As with robust optimization, parameters on security-level and probability measures are to be given to the model, which is a difficult task in practice.

We contribute to the existing literature by considering inflow uncertainty in the short-term hydropower model. Few papers have looked specifically in stochastic short-term models and we extend the modeling proposed in [91] to consider uncertain inflows. For the producer, it is interesting to consider a stochastic model since it gives a production plan for the whole planning horizon. Applying the theory outlined in [82], we also detail/provide a nonparametric scenario generation approach that relies on the information in the history of inflows. We expand [91] by introducing stochasticity to both the loading and unit commitment problems.

**The paper is organized as follows** Section 5.2 presents data available for inflows. Section 5.3 describes the method to generate scenario trees. Section 5.4 gives an overview of the short-term hydropower problem and details the optimization models. Numerical results are presented in Section 5.5 and final remarks are presented in Section 5.6.

## 5.2 Scenario fan of inflows

This section presents the data available for the inflows. In the province of Quebec in Canada, consumers and producers of hydroelectric energy, except Hydro-Quebec, are not allowed to bid on the spot markets [22]. The province-owned integrated utility performs all power market activities. Hence, only uncertainty related to inflows in the reservoirs is considered in this paper.
Before presenting the method for generating the scenario trees used in the optimization models, we describe the available data sets. Precipitation forecasts are obtained from Environment Canada [1]. A 7 day deterministic precipitation forecast is issued. The 7 day forecast is split in two groups : the first 3 and the last 4 days. We make the assumption that the error for both groups is independent from a meteorological point of view, as the correlation in precipitations between days is negligible. This assumption is motivated by the high variability of the weather conditions on our watershed from one day to the next. For example, we could have a few days of snow, followed by no precipitations then a few days of rain. The last 15 years of historical data of precipitation forecasts is searched for a given number (a) of precipitation forecasts that are the closest, in precipitation forecast (mm) to the first 3 days, and they are retained. The same is conducted for the second group. Since the error is assumed independent, the scenarios found for the first and second group are mixed and matched to create a² precipitation scenarios for the first 7 days. Note that the actual realizations of precipitation on these days are used as scenarios. Then, considering that the forecast has no value after 7 days, the 62 years of available history of realizations is appended to all of the scenarios for the first 7 days with a = 7, yielding a total of a² × 62 = 3038 scenarios of precipitation for 30 days of prevision. Then, these precipitation scenarios are given as input to the CEQUEAU hydrological model [25] which outputs inflow previsions for the reservoirs.

Figure 5.1 illustrates this process. The goal of the scenario tree generation method, in Section 5.3, is to create a scenario tree from the scenario fan of inflows.

**5.3 Scenario tree generation**

The method chosen to construct a scenario tree suitable for the stochastic optimization is taken from [82, 81]. The method is applied to real hydropower data. First, the structure of the scenario tree is fixed, then stochastic approximation is used to improve the values of inflow of the nodes, considering all the data available for every approximation. Improvement goes on until a convergence criterion, based on the nested distance and explained in Section 5.3.4, is reached.

**5.3.1 Fixing the initial scenario tree structure : k-means clustering**

The number of stages and the number of nodes per stage of the tree are fixed initially. Aggregation is necessary since the scenario tree structure can be different from the data available. The aggregation is straightforward : values of inflows for each day are summed up. K-means clustering [68] is used to partition the data paths into clusters in order to assign
<figure><img src="image_18.png" /><figcaption>Figure 5.1 Building inflow scenarios from a 7 day deterministic precipitation forecast.</figcaption></figure>

initial values to the scenario tree nodes. Note that initially no probabilities are allocated to the nodes : simply values for the nodes. This clustering method minimizes the distance from every data point to the mean of the cluster to which it belongs. As an example, the k-means algorithm is applied to the 3038 inflow scenarios to form a scenario tree which has a structure as shown in Figure 5.3b.

### 5.3.2 Improvement of the clusters

The method to improve the scenario tree nodes consists of two steps. First, from the initial data paths, a random data path, that is not in the paths available, is generated using density estimation. Next, the distance between this random path and the closest state of the scenario tree nodes is minimized in a stochastic approximation step in order to improve the tree. This method is repeated for a given number of iterations and is explained in what follows.

#### Step I : density estimation

In order to generate a new random path, kernel density estimation is used. We generate a random path that is close to the distribution of the data paths and conditional on previous stages. To do so, the conditional probability density function is estimated. For each stage of the desired scenario tree structure, a value of inflow is generated that is close in distribution to all of the data paths and incidental to the past.

A random path $\xi_k^d = (\xi_1^d, \dots, \xi_K^d)^T$ is to be generated using available data paths $X_{ik}^d = (X_{i1}^d, \dots, X_{iK}^d)^T$ where $i$ is the index of available data paths, $d$ is the dimension and $K$ is the
number of stages. The conditional density estimator [81] is :

$$ \hat{f}(\xi_k | \xi_1, \dots, \xi_{k-1}) = \sum_{i=1}^{n} \prod_{j=1}^{k-1} \frac{\kappa\left(\frac{\xi_j - X_{ij}}{h_j}\right)}{\sum_{m=1}^{n} \kappa\left(\frac{\xi_j - X_{jm}}{h_j}\right)} \times \kappa\left(\frac{\xi_k - X_{ik}}{h_k}\right) \times \frac{1}{h_k}, \quad (5.1) $$

where the dimension $d$ is dropped for clarity, $n$ is the number of available data paths, $\kappa$ is the kernel and $h_j$ and $h_k$ are the bandwidths.

The analytical representation of the actual distribution is not computed, as only samples from Equation (5.1) are necessary which can be generated quickly. In practice, this is achieved by assigning weights to every data path available. The closer the observation is to the path, the higher is the weight. For every stage from 1, ..., $k-1$, the weights of the data path at each stage are multiplied. With these weights calculated, a value of inflow is to be generated at stage $k$.

<figure><img src="image_19.png" /><figcaption>Figure 5.2 Generation of a random path based on three available data paths of inflows. The generated value of inflow for stage 1 is shown with a star marker.</figcaption></figure>

To illustrate refer to Figure 5.2. There are three data paths of inflow. The random value of inflow has been generated for stage 1 and is located with a star marker. From there, a value of inflow is to be generated for subsequent stages, always conditional on the past. As per the figure, it is necessary to find a value of inflow at stage 2 that is consistent with the conditional distribution. Therefore, weights are calculated as follows, in this case for stage $k$:

$$ w_i(\xi_1, \dots, \xi_{k-1}) = \prod_{j=1}^{k-1} \frac{\kappa\left(\frac{\xi_j - X_{ij}}{h_j}\right)}{\sum_{m=1}^{n} \kappa\left(\frac{\xi_j - X_{jm}}{h_j}\right)}, \quad (5.2) $$
where $\sum_{i=1}^{n} w_i = 1$ and $w \ge 0$.

The value of inflow $\xi_k$ at stage $k$ is generated as follows. A data path with index $i^*$ is chosen randomly among the available data paths at stage $k-1$ to satisfy

$$ \sum_{i=1}^{i^*-1} w_i(\xi_1, \dots, \xi_{k-1}) \le \mathrm{rand}_u \le \sum_{i=1}^{i^*} w_i(\xi_1, \dots, \xi_{k-1}), \quad (5.3) $$

where $\mathrm{rand}_u$ is chosen from the uniform random distribution on the interval $[0, 1]$. The cumulative sum of the weights leads to a high probability of picking a data path near an observation.

The value of inflow $\xi_k$ is obtained by setting the value at stage $k$

$$ \xi_k = X_{i^*k} + \mathrm{rand}_{\kappa_{h_k}}, \quad (5.4) $$

where $\mathrm{rand}_{\kappa_{h_k}}$ is a random value sampled from the kernel estimator using the composition method [82].

This newly generated inflow value is according to the distribution of density of the current stage and dependent on the history of all the data paths.

Referring again to Figure 5.2, weights are calculated for the 3 data paths as per Equation (5.2). Then, a data path is chosen randomly at stage 1 and the solid line has a high probability of being picked. Consider it is the case. To generate the value of inflow at stage 2, the value of the solid line at stage 2 is perturbed randomly. This method is then repeated at each stage in order to generate a random data path and is represented on Figure 5.3a with a thick dashed line.

It is shown that the choice of the kernel does not have an important effect on the density estimation [58]. Hence, in this paper, the logistic kernel is used :

$$ \kappa(\xi) = \frac{1}{e^{\xi} + 2 + e^{-\xi}}. \quad (5.5) $$

The bandwidth is the smoothing factor applied to the estimation of the density. Silverman's rule of thumb [94] is employed to determine the optimal bandwidth :

$$ h_k = \sigma(X_{ik})n^{-\frac{1}{d+4}} = \sigma(X_{ik})n^{-\frac{1}{7}}, \quad (5.6) $$

where $n$ is the number of data paths, $d$ is the dimension and $\sigma$ is the standard deviation. In this paper, $d=3$ because there are three values of inflows per scenario tree node, representing
three different reservoirs.

## Step II : stochastic approximation

Once the new random path of inflows is generated, a stochastic approximation step is conducted. This step allows to update the value of some scenario tree states. During this step, a scenario from the scenario tree, more precisely a path of nodes in the scenario tree is identified. This path of nodes in the scenario tree minimizes the Wasserstein distance $W$ between the randomly generated path during Step I of the algorithm, found in Section 5.3.2, and current scenario tree nodes values.

The Wasserstein distance is minimized as follows :

$$W^2 = \min_{\omega \in \Omega} \sum_{k=1}^{K} ||\Gamma(\omega) - \xi_k||^2, \qquad (5.7)$$

where $\Omega$ are the scenario tree paths, $\Gamma(\omega)$ are the states corresponding to the nodes in the path $\omega$ in the scenario tree, from the set of all possible scenarios $\Omega$, and $\xi_k$ is the value of inflow generated randomly at stage $k$. Referring to Figure 5.3b, $\Omega = \{(1, 2, 3, 5), (1, 2, 3, 6), (1, 2, 4, 7), (1, 2, 4, 8)\}$. Equation (5.7) allows to find this path of nodes and is identified as nodes (1, 2, 4, 8) on Figure 5.3b.

To achieve this, a stochastic gradient descent method that minimizes the nested distance is used. Starting from the root of the scenario tree, $W$ is computed for the children node. The children node with the smallest value of $W$ becomes the parent node. $W$ is then computed for the children node of the new parent node and so on until a leaf node has been reached.

The identified path of scenario tree nodes values $\Gamma(\omega)$ that minimizes the Wasserstein distance for the current stochastic approximation iteration $p = 1, 2, \dots$ is updated in the following manner :

$$\Gamma(\omega)_{p+1} = \Gamma(\omega)_p - \alpha_p \nabla W_p, \qquad (5.8)$$

where $\Gamma(\omega)$ are the values of the scenario tree nodes to improve, $\alpha_p$ is the step-size and $\nabla W_p$ the gradient of the distance.

The step-size $\alpha_p = \frac{1}{(p+30)^{3/4}}$, where $p$ is the stochastic approximation iteration, is chosen since it is shown that the method will converge [81] given $\alpha_p > 0$, $\sum_p \alpha_p = \infty$ and $\sum_p (\alpha_p)^2 < \infty$.

As an illustration, consider one iteration of the algorithm and refer to Figure 5.3. First, a random data path of inflows is generated using kernel density estimation. This can be seen on Figure 5.3a : it is the thick dashed line. The Wasserstein distance between this new generated path of inflows and the current values of the scenario tree nodes is minimized and
<figure><img src="image_20.png" /><figcaption>Figure 5.3 Illustration of the 2 steps of the algorithm. Generation of a random path of inflows from available data paths of inflows and stochastic approximation to improve the value of some scenario tree nodes.</figcaption></figure>

a path of nodes in the scenario tree is retrieved for potential improvement. The path of nodes minimizing this distance is shown on Figure 5.3b. Hence, the value of the inflows for the thick nodes, which are 1, 2, 4 and 8 will be improved using Equation (5.8).

### 5.3.3 Probabilities

During the first stochastic approximation iteration, assigned probabilities of the nodes are 0 since, as explained in Section 5.3.1, the scenario tree is initialized with values for the nodes only.

Node probabilities are updated at every stochastic approximation iteration. A counter is assigned to each node and initialized at 0. Every time a path of nodes minimizing the Wasserstein distance is retrieved, the corresponding counters of the nodes in this path are incremented by 1.

Once the stochastic approximation iterations are completed, probabilities are computed by dividing the counter value by the number of stochastic approximation iterations, which yields sums of child nodes probabilities equal to 1, as in Figure 5.4.

In a multistage stochastic program, each path from the root to a leaf node represents a scenario. The unconditional probabilities of a scenario is obtained by multiplying the unconditional probabilities of all the nodes in the scenario, yielding probabilities $\pi_j$, where $j$ is the scenario in Figure 5.4.

An interesting feature of the scenario tree generation method is that the extreme (low and
<figure><img src="image_21.png" /><figcaption>Figure 5.4 A scenario tree with node probabilities (over the node) and scenario *i* probabilities (indicated with $\pi_i$)</figcaption></figure>

high) scenarios are accounted for, according to their occurrence in the historical data set. The law of large numbers insures that the probabilities are asymptotically consistent.

### 5.3.4 Termination criterion

The scenario tree generation algorithm terminates when the nested distance has converged to a certain $\epsilon$ for the 10 last iterations. Thus, Step I and Step II of the algorithm are repeated until convergence is obtained. Depending on the inflow forecasts, the number of iterations to converge varies. As an example, for a given test case, it took an average of 1038 iterations for the method to converge and generate one scenario tree.

The main advantage of the scenario tree generation method presented in this section is that all of the data paths are used at every iteration to improve the value of the scenario tree nodes. By doing so, the underlying discrete distribution of the available data paths, approximated by a scenario tree, is improved consistently with the data. The scenario trees are prepared before the optimization is conducted.
## 5.4 Stochastic short-term hydropower model

The two-phase deterministic optimization models taken from [91] are updated to consider stochastic inflows. This section presents the modeling of the short-term problem as well as the mathematical formulations.

### 5.4.1 Modeling of the short-term problem

The modeling of the problem considers head-dependency, as well as efficiencies of each turbine. Power $P(kW)$ produced by a single turbine is defined as

$$P(h_n, Q) = \eta(Q) \times G \times Q \times h_n(Q_{tot}, v), \qquad (5.9)$$

where $G$ is the gravitational acceleration ($m/s^2$), $Q$ is the unit water flow and $Q_{tot}$ is the total water flow ($m^3/s$), $\eta(Q)$ is the efficiency of the turbine and $h_n$ is the net water-head ($m$). The net water-head is a function of the forebay elevation $h_f(m)$, the tailrace elevation $h_t(m)$ and losses in the penstock $\varphi(m)$ that is given by :

$$h_n(Q_{tot}, v) = h_f(v) - h_t(Q_{tot}) - \varphi(Q_{tot}), \qquad (5.10)$$

where $v$ is the volume of the reservoir ($hm^3$). For notational purposes and since there is a relation between net water head and volume, we consider that power is a function of the volume and water flow. We propose a modeling with combinations of units instead of single units. To achieve this, a dynamic programming algorithm, where each sub-problem is a turbine, is used to calculate the power produced by a combination of units. As an example, if a power plant has a total of 5 turbines and requires three active turbines, there is a total of 10 combinations of 3 turbines, 5 combinations of 4 turbines and 1 combination of 5 turbines. Water flows are discretized and the dynamic programming algorithm is executed for each possible combinations, 16 in this case, for each power plant and discretizations of reservoir volumes and water flows.

#### Dynamic programming algorithm

The objective of the problem is to maximize the power output and it is found recursively. Given state $s^j$, the dynamic programming algorithm seeks to choose decision variables $q^j$ that solves :

$$f^{*j}(s^j) = \max_{q^j} P(s^j, v) + f^{*j+1}(s^j - q^j), \qquad (5.11)$$
where $j = n - 1, n - 2, \dots, 1$, $n$ is the number of turbines at the power plant, $s^j \in \{1, 2, \dots, r\}$ is the remaining water to dispatch given the number of discretizations $r$ and $q^j \in \{1, 2, \dots, \min\{\bar{q}^j, Q\}\}$ the water flow with $\bar{q}^j$ maximum water flow. The optimal water flow is $q^{*j} = s^j$ that maximizes $f^{*j}(s^j)$. For $j = n$, the optimal power output is given by $f^{*j}(s^j) = P(s^j)$.

## Maximum power output surfaces

We then build a set of surfaces of the maximum power output for each power plant. For a plant with 5 turbines requiring at least 3 working, three surfaces are built, more precisely one for 3 turbines working, one for 4 turbines working and one for 5 turbines working. The maximum power output for every possible combination of number of working turbines is retained for every discretization of volume and water flow. Such surfaces can be viewed in Figure 5.5. To

<figure><img src="image_22.png" /><figcaption>Figure 5.5 Maximum output surfaces</figcaption></figure>

obtain them, the dynamic programming algorithm is executed for every number of turbines in the combination, every discretization of the reservoir volume, every discretization of the water flow and every power plant. The surfaces of maximum power outputs are then modeled using polynomial equations in the objective-function of the optimization problem. Modeling of the hydropower production functions is done by constraining these functions with two surfaces.

One constrains the water discharge without spill and the other constrains the water discharge with spill. Therefore, in the optimization problem, we have only one variable of the water discharge, which combines processed and spilled. But since we are constraining with the two
different surfaces and that we are maximizing, the model will try and avoid spill since it
reduces power production. When redistributing the total water discharge to the plant, we
can then see if there is some spill or not. There is an upper bound on the water discharge,
which is the maximum including spill.

A two-phase optimization strategy is used to penalize the startup of turbines. The first phase, namely the loading problem, optimizes values of water discharges, volumes and number of turbines in the combination for every plant and node. The second phase, namely the unit commitment problem, uses the solution of the first optimization model to determine the exact combination of turbines working at each plant and node in the scenario tree. Startups of turbines are penalized with a fixed cost. Multistage stochastic models are developed for both optimization phases, in order to consider uncertainty in the inflows of the reservoirs.

5.4.2 Phase I : loading problem

Optimization variables of this nonlinear stochastic multistage mixed integer problem are water flows, volumes and number of working turbines, for each node and plant in the scenario tree. There is only one variable for the water flow, but it includes processed and spilled water. We have shown [91] that relaxing the variables that determine the number of working turbines leads to an integer solution. Therefore, we solve a nonlinear stochastic multistage continuous problem, as the coefficients of the matrix of constraints are also totally unimodular given the stochastic version of the model.

The objective is to maximize total energy production in stage 0, expected energy production in future stages and expected value of the water remaining in the reservoir at the end of the planning horizon :

$$
\max_{y,q,v} \sum_{c \in C} \sum_{s=1}^{n_0^c} \chi_{s0}^c y_{s0}^c \zeta_0 + \sum_{c \in C} \sum_{j \in K} \pi_j^c \left( \sum_{i \in N_j} \sum_{s=1}^{n_i^c} \chi_{si}^c y_{si}^c \zeta_i + \sum_{p \in E_j} \Phi_p^c (v_p^c) \right) \quad (5.12)
$$

$$
\begin{array}{ll}
\text{subject to :} & \chi_{si}^c \le \Psi_s^{Ac}(q_i^c, v_i^c), \quad \forall c \in C, \forall i \in N, \\
& \qquad \forall s \in \{1, 2, \dots, n_i^c\}, \hfill (5.13)
\end{array}
$$

$$
\begin{array}{ll}
\chi_{si}^c \le \Psi_s^{Bc}(q_i^c, v_i^c), & \forall c \in C, \forall i \in N, \\
& \forall s \in \{1, 2, \dots, n_i^c\}, \tag{5.14}
\end{array}
$$

$$
\begin{equation}
\begin{split}
\delta_i^c = v_{i+1}^c & - v_i^c + \gamma w_i q_i^c \\
& - \sum_{m=1}^{u^c} \gamma w_m q_i^m,
\end{split}
\tag{5.15}
\end{equation}
$$
$$ \sum_{s=1}^{n_i^c} y_{si}^c \le 1, \quad \forall i \in N, \forall c \in C, \qquad (5.16) $$

$$ y_{s0}^c = \hat{y}_{s0}^c, \quad \forall s \in \{1, 2, \dots, n_i^c\}, \forall c \in C, \forall i \in N, \qquad (5.17) $$

$$ v_{min}^c \le v_i^c \le v_{max}^c, \quad \forall i \in N, \forall c \in C, \qquad (5.18) $$

$$ q_{min}^c \le q_i^c \le q_{max}^c, \quad \forall i \in N, \forall c \in C, \qquad (5.19) $$

$$ q_i^c \ge 0, \quad \forall i \in N, \forall c \in C, \qquad (5.20) $$

$$ v_i^c \ge 0, \quad \forall i \in N, \forall c \in C, \qquad (5.21) $$

$$ y_{si}^c \ge 0, \quad \forall s \in \{1, 2, \dots, n_i^c\}, \forall i \in N, \forall c \in C. \qquad (5.22) $$

Hydropower production surfaces are constrained by (5.13)-(5.14). Water balance constraints are represented by (5.15) and the choice of a single number of active turbines is shown in (5.16). Constraints (5.17) are the initial number of active turbines while constraints (5.18)-(5.19) are the bounds on reservoir volumes and water discharges. Finally, constraints (5.20)-(5.22) impose nonnegativity. Index $i+1$ of variables $v_{i+1}^c$ takes the value of the node in the set $N_j$. For example, if $N_j = \{1, 3, 5, 7\}$ and $i = 1, i+1 = 3$ since the index takes the value of the node at position $i=2$ in the set $N_j$.

The above short-term loading problem is described in more details in [91]. We now show how to integrate a water-value function for the remaining water at the end of the planning horizon.

**Water-value function** The water-value function is the expected energy production in the future at the end of the planning horizon. In a deterministic framework, inflows are known with certainty, thus volume in the reservoir at the end of the horizon is easier to determine. In a stochastic framework, it is not possible to give a goal for the volume at the end of the horizon since it may not be feasible for every scenario. On the other hand, neglecting this feature will cause the optimization to empty the reservoir at the end of the horizon, since the objective is to maximize energy. Hence, maximizing the expected value of future energy production, or water-value function, will prevent the optimization of doing this. The water-value functions are computed with a stochastic dynamic algorithm [33] at Rio Tinto. A planning horizon of one year, with weekly time steps is used.

### 5.4.3 Phase II : unit commitment

This linear stochastic multistage integer model is solved using solution found in Phase I. The purpose of this model is to determine the on-off schedule of the turbine combinations (found in Phase I). Given water flows and reservoir volumes found in the loading problem, the dynamic programming algorithm is used to calculate power outputs for every possible
combination of turbines, given the number of working turbines found in Phase I, and are stored in parameter $\beta_{li}^c$. The model maximizes the energy production and penalizes turbine startups. Initial combination of turbines working at stage 0 is given in $\hat{x}_{l0}^c$.

The objective is to maximize energy production at stage 0 and future energy production and penalize startup of turbines at stage 0 as well as future expected startups:

$$ \max_{x,d} \sum_{c \in C} \zeta_0 \left( \sum_{l=1}^{n_0^c} \beta_{l0}^c x_{l0}^c - \sum_{t=1}^{T^c} d_{l0}^c \theta \right) + \sum_{j \in K} \sum_{c \in C} \left( \pi_j^c \left( \sum_{i \in N_j} \zeta_i \left( \sum_{l=1}^{n_i^c} \beta_{li}^c x_{li}^c - \sum_{t=1}^{T^c} d_{li}^c \theta \right) \right) \right) \quad (5.23) $$

$$ \text{subject to :} \quad \sum_{l=1}^{n_i^c} x_{li}^c = 1, \quad \forall i \in N, \forall c \in C, \quad (5.24) $$

$$ x_{li}^c f_{lit}^c - x_{li-1}^c f_{li-1t}^c \le d_{ti}^c, \quad \forall l \in \{1, 2, \dots, n_i^c\}, \forall i \in N_j, \forall j \in K, \forall c \in C, \\ \forall t \in \{1, 2, \dots, T^c\}, \quad (5.25) $$

$$ x_{l0}^c = \hat{x}_{l0}^c, \quad \forall l \in \{1, 2, \dots, n_i^c\}, \forall i \in N, \forall c \in C, \quad (5.26) $$

$$ x_{li}^c, d_{it}^c \in \mathbb{B}, \quad \forall l \in \{1, 2, \dots, n_i^c\}, \forall i \in N, \forall t \in \{1, 2, \dots, T^c\}, \forall c \in C. \quad (5.27) $$

The choice of a single turbine combination is given by (5.24). Constraints that allow to penalize a startup by flagging them is shown in constraints (5.25). The initial combinations are given in (5.26) and imposition of binary variables are constraints (5.27). Index $i-1$ of parameters $f_{li-1t}^c$ takes the value of the node in the set $N_j$. For example, if $N_j = \{1, 3, 5, 7\}$ and $i=4, i-1=5$ since the index takes the value of the node at position $i=3$ in the set $N_j$.

This two-phase optimization process allows to find a solution efficiently. Also, even though an approximation of the energy produced is conducted in the first phase, the actual energy production is retrieved in the second phase, seeing that the actual hydropower production functions are used to compute the actual energy production given a water discharge and volume, which are solutions of the first phase.

## 5.5 Results

This section details the system on which the stochastic hydropower models are tested and results are presented.

### 5.5.1 Hydroelectric system studied

The hydroelectric system studied is situated in the Saguenay Lac-St-Jean region in the province of Quebec, Canada and is owned by Rio Tinto. For the purpose of this paper, three
hydroelectric plants, which are Chute-du-Diable, Chute-Savane and Isle-Maligne are considered. The two first plants have 5 turbines each and the latter has 12. Figure 5.6 represents the system studied. Triangles represent reservoirs and squares power plants.

<figure><img src="image_23.png" /><figcaption>Figure 5.6 Hydroelectric system studied.</figcaption></figure>

Chute-du-Diable, Chute-Savane and Isle-Maligne plant reservoirs are quite small, respectively 452 $hm^3$, 119 $hm^3$ and 171 $hm^3$. In the optimization model, there is no water value function associated to these plants since they have small reservoirs. Instead, a full reservoir constraint at the last period is imposed as a goal in the model. The only water-value function used is for the Lac-St-Jean reservoir, therefore volume of this reservoir at the last period is an optimization variable. The capacity of this reservoir is of 5596 $hm^3$. Water flow in *Petite décharge* is limited by a function dependent on the volume of Lac-St-Jean.

### 5.5.2 Rolling-horizon procedure

A rolling-horizon methodology is retained to validate the optimization models developed in this paper. The planning horizon of the rolling-horizon is of 31 days. For every day of the rolling-horizon, the forecast is for 30 days. For day 1 of the rolling-horizon, previsions are from days 2 to 31, for day 2 of the rolling-horizon, previsions are from days 3 to 32, and so on. The stochastic optimization models presented in section 5.4 are solved every day, but only the solution for the first-stage is retained. Forecasts are updated daily. Once the forecast is updated, the scenario tree is generated for the corresponding day. The two-phase optimization process is launched and the first-stage solution is retained, that is : volume, water discharge and turbine combination. Then, considering the actual realization of the inflow, the water balance constraints are used to determine the actual volume of the reservoirs at the end of the period. More precisely, the water discharge from the optimization is combined with the
actual realization of the inflow in order to calculate the reservoir volumes. The same process is
repeated for the 31 days. In the end, a production plan for 31 days is available, which consists
of the reservoir volumes, total water discharges at the plants and turbine combinations in
use. See [3] for a different approach to rolling-horizon evaluation of short-term hydropower
operation.

The solution obtained from the scenario tree generation is compared to the solution obtained
from the median scenario of the inflows. Therefore, we compare our method to a rolling
median. Every day, the median scenario is found throughout all available scenarios and a
scenario tree of 1 node per stage is solved in a deterministic fashion.

5.5.3 Numerical results

The scenario tree generation method is coded in Matlab [72]. The optimization models are coded using AMPL [48]. The optimization software for the loading problem, which is the relaxation of a nonlinear mixed-integer problem, is IPOPT [105] and the unit commitment model, a linear integer problem, is solved with XPRESS [2].

Six test cases, which consist of monthly periods are available. The biggest problems to
solve have 7 stages with 48 scenarios, 1123 nonlinear variables, 33 linear variables and 1237
constraints for the loading problem and 3475 binary variables and 825 constraints for the
unit commitment problem.

Different stages, more precisely 5, 6 or 7, as well as different number of scenarios, namely 16, 32 or 48 are tested.

Computational time

The average time to construct the scenario tree and to optimize is shown on Figure 5.7. The average time is in seconds, for a single day in the rolling-horizon procedure, more precisely for one problem including construction of the scenario tree and optimization of the two-phase process. It takes less than 5 seconds to build the scenario trees for all test cases, while the optimization requires more time given higher numbers of scenarios. Less than 42 seconds, for a single day in the rolling-horizon are necessary to construct the scenario tree and optimize the two-phase process, which is acceptable in the real operating environment.

The current implementation of the scenario tree generation method and optimization is
tested on three cascaded hydropower plants. For this specific producer, the whole hydropower
systems consists of five hydropower plants, therefore calculation time would be acceptable
for the whole system. Considering another system of, for example, 50 hydropower plants,
<figure><img src="image_24.png" /><figcaption>Figure 5.7 Average computational time of scenario tree generation and optimization for one day in the rolling-horizon.</figcaption></figure>

the actual method would take approximately 350 minutes. The proposed method in this paper is applicable to a larger system, probably by decomposing the system in smaller subsystems. To do so, the system is to be studied and depending on its configuration, distances between plants and others, modeled in an acceptable manner. Depending on the scope of the application, the calculation time may or may not be satisfactory. If a producer does not mind solving a 7 hour model every day, then the computational time is satisfactory. In order to diminish computing time, an avenue is to solve the model for a given number of days then weeks. In this way, the number of variables is greatly reduced and so is the computing time. This model is applicable to a larger hydropower system, but it would be necessary to decompose the system in sub-systems and review the modeling to diminish the number of optimization variables, given a producer requiring fast computational time.

## Results

Table 5.1 illustrates the difference in energy, in TWh, produced throughout the 31 days rolling-horizon combined with the value of water remaining in the reservoir at the end of the planning horizon. This implies that the difference in energy can be compared to annual production but absolute numbers are unfortunately not thus interpretable. A positive value indicates the scenario tree method produces more than the median scenario and a negative value indicates the contrary. For 4 of the test cases, the stochastic solution produces more energy. For 1 test case, the median scenario solution produces more energy. Finally, for the August case, the stochastic solution produces more energy with a 5 stage or 6 stage scenario
<table><caption>Table 5.1 Results for 6 test cases (5 are data sets from the year 2011 and 1 from 2010). Energy produced by the stochastic solution and the median scenario rolling-horizon is given. Also, the difference in energy between both solutions is shown.</caption><thead><tr><th rowspan="2">Nb. Sc.</th><th colspan="3">June 2011</th><th colspan="3">July 2011</th><th colspan="3">August 2011</th></tr><tr><th>Stoch. (TWh)</th><th>Median (TWh)</th><th>Diff. (TWh)</th><th>Stoch. (TWh)</th><th>Median (TWh)</th><th>Diff. (TWh)</th><th>Stoch. (TWh)</th><th>Median (TWh)</th><th>Diff. (TWh)</th></tr></thead><tbody><tr><td colspan="10" style="text-align:center; font-weight:bold;">5 stages</td></tr><tr><td>16</td><td>804.5143</td><td>804.0265</td><td>0.4878</td><td>740.2678</td><td>740.0631</td><td>0.2047</td><td>710.1115</td><td>710.0795</td><td>0.0320</td></tr><tr><td>32</td><td>804.7050</td><td>804.0251</td><td>0.6799</td><td>740.2783</td><td>740.0631</td><td>0.2152</td><td>710.1108</td><td>710.0794</td><td>0.0314</td></tr><tr><td>48</td><td>804.6894</td><td>804.0249</td><td>0.6645</td><td>740.2496</td><td>740.0631</td><td>0.1865</td><td>710.0988</td><td>710.0794</td><td>0.0194</td></tr><tr><td colspan="10" style="text-align:center; font-weight:bold;">6 stages</td></tr><tr><td>16</td><td>804.5059</td><td>804.1495</td><td>0.3564</td><td>740.2698</td><td>740.0665</td><td>0.2033</td><td>710.0783</td><td>710.0733</td><td>0.0050</td></tr><tr><td>32</td><td>804.6796</td><td>804.1479</td><td>0.5317</td><td>740.2652</td><td>740.0665</td><td>0.1987</td><td>710.1139</td><td>710.0733</td><td>0.0406</td></tr><tr><td>48</td><td>804.6715</td><td>804.1481</td><td>0.5234</td><td>740.2608</td><td>740.0665</td><td>0.1943</td><td>709.9826</td><td>710.0732</td><td><strong>-0.0906</strong></td></tr><tr><td colspan="10" style="text-align:center; font-weight:bold;">7 stages</td></tr><tr><td>16</td><td>804.5171</td><td>804.0881</td><td>0.4290</td><td>740.2676</td><td>740.0578</td><td>0.2098</td><td>710.0693</td><td>710.0867</td><td><strong>-0.0174</strong></td></tr><tr><td>32</td><td>804.7166</td><td>804.0881</td><td>0.6285</td><td>740.2566</td><td>740.0578</td><td>0.1988</td><td>710.0732</td><td>710.0867</td><td><strong>-0.0135</strong></td></tr><tr><td>48</td><td>804.7063</td><td>804.0879</td><td>0.6184</td><td>740.2686</td><td>740.0578</td><td>0.2108</td><td>710.0806</td><td>710.0867</td><td><strong>-0.0061</strong></td></tr><tr><th rowspan="2">Nb. Sc.</th><th colspan="3">September 2010</th><th colspan="3">September 2011</th><th colspan="3">October 2011</th></tr><tr><th>Stoch. (TWh)</th><th>Median (TWh)</th><th>Diff. (TWh)</th><th>Stoch. (TWh)</th><th>Median (TWh)</th><th>Diff. (TWh)</th><th>Stoch. (TWh)</th><th>Median (TWh)</th><th>Diff. (TWh)</th></tr><tr><td colspan="10" style="text-align:center; font-weight:bold;">5 stages</td></tr><tr><td>16</td><td>729.5792</td><td>729.3811</td><td>0.1981</td><td>733.0375</td><td>731.6799</td><td>1.3576</td><td>704.7842</td><td>704.8494</td><td><strong>-0.0652</strong></td></tr><tr><td>32</td><td>729.5841</td><td>729.3821</td><td>0.2020</td><td>733.0530</td><td>731.6799</td><td>1.3731</td><td>704.7847</td><td>704.8494</td><td><strong>-0.0647</strong></td></tr><tr><td>48</td><td>729.5810</td><td>729.3804</td><td>0.2006</td><td>733.0818</td><td>731.6799</td><td>1.4019</td><td>704.7877</td><td>704.8496</td><td><strong>-0.0619</strong></td></tr><tr><td colspan="10" style="text-align:center; font-weight:bold;">6 stages</td></tr><tr><td>16</td><td>729.5856</td><td>729.3917</td><td>0.1939</td><td>732.9971</td><td>731.7773</td><td>1.2198</td><td>704.7690</td><td>704.8636</td><td><strong>-0.0946</strong></td></tr><tr><td>32</td><td>729.5779</td><td>729.3929</td><td>0.1850</td><td>733.0188</td><td>731.7773</td><td>1.2415</td><td>704.7928</td><td>704.8636</td><td><strong>-0.0708</strong></td></tr><tr><td>48</td><td>729.5800</td><td>729.3924</td><td>0.1876</td><td>733.0937</td><td>731.7774</td><td>1.3163</td><td>704.7326</td><td>704.8634</td><td><strong>-0.1308</strong></td></tr><tr><td colspan="10" style="text-align:center; font-weight:bold;">7 stages</td></tr><tr><td>16</td><td>729.5854</td><td>729.4151</td><td>0.1703</td><td>732.9428</td><td>731.9647</td><td>0.9781</td><td>704.7608</td><td>704.8566</td><td><strong>-0.0958</strong></td></tr><tr><td>32</td><td>729.5775</td><td>729.4156</td><td>0.1619</td><td>732.9599</td><td>731.9647</td><td>0.9952</td><td>704.7879</td><td>704.8566</td><td><strong>-0.0687</strong></td></tr><tr><td>48</td><td>729.5834</td><td>729.4139</td><td>0.1695</td><td>732.9702</td><td>731.9648</td><td>1.0054</td><td>704.7873</td><td>704.8567</td><td><strong>-0.0694</strong></td></tr></tbody></table>

tree, and the median scenario with a 7 stage. For the 4 test cases for which the scenario tree produces more energy than the median scenario, average improvements are 0.0679812% for June, 0.0273551% for July, 0.1620522% for September 2011 and 0.0251653% for September 2010. Despite the low percentages, this represents significant savings for the producer. As an example, the current value of a 1 GWh improvement, in the province of Quebec, is around 20,000$. Therefore, for June, the 0.0679812% higher production translates into 10,932,489$.
## In-sample stability test

An in-sample stability test allows to verify if the scenario tree generation method is consistent. It is taken from [62]. Since the scenario tree is generated from random samples, one wants to verify if the solution given by the optimization, with a different scenario tree each time, give more or less the same solution. If so, then the scenario tree method is consistent.

As an example, July 2011 and June 2011 data sets were chosen for this verification. For both data sets, 6 scenario trees were generated with the same number of stages and scenarios. Then, the optimization was conducted on all of these scenario trees to verify the effect on the objective function value. Table 5.2 gives, for these two data sets and 6 instances each, the values of the objective function, for the scenario tree and median scenario methods.

<table><caption>Table 5.2 Objective function values for 6 random scenario trees with the same number of stages and scenarios, on two data sets.</caption><thead><tr><th>Data</th><th>Inst.</th><th>Stoch.<br>(TWh)</th><th>Median<br>(TWh)</th><th>Diff.<br>(TWh)</th></tr></thead><tbody><tr><td rowspan="6">July</td><td>1</td><td>740.2652</td><td>740.0665</td><td>0.1987</td></tr><tr><td>2</td><td>740.2759</td><td>740.0665</td><td>0.2094</td></tr><tr><td>3</td><td>740.2725</td><td>740.0665</td><td>0.2060</td></tr><tr><td>4</td><td>740.2581</td><td>740.0665</td><td>0.1916</td></tr><tr><td>5</td><td>740.2799</td><td>740.0665</td><td>0.2134</td></tr><tr><td>6</td><td>740.2878</td><td>740.0665</td><td>0.2213</td></tr><tr><td rowspan="6">June</td><td>1</td><td>804.6715</td><td>804.1481</td><td>0.5234</td></tr><tr><td>2</td><td>804.6707</td><td>804.1484</td><td>0.5223</td></tr><tr><td>3</td><td>804.6709</td><td>804.1474</td><td>0.5235</td></tr><tr><td>4</td><td>804.6824</td><td>804.1489</td><td>0.5335</td></tr><tr><td>5</td><td>804.6769</td><td>804.1486</td><td>0.5283</td></tr><tr><td>6</td><td>804.6571</td><td>804.1472</td><td>0.5099</td></tr></tbody></table>

Results show that the scenario tree generation method is consistent, as the difference between the objective functions of the stochastic and median scenario methods present slight variations. For the July test case, the median is 0.2077 TWh, the mean 0.2067 TWh and the variance 0.9308 TWh and for the June test case, the median and the mean are 0.5235 TWh and the variance 0.0516 TWh.

## Interpretation of the results

The following figures illustrate the 31 day rolling-horizon backtesting solution more precisely: water discharge and reservoir levels for the power plants and reservoirs studied in this paper.
<figure><img src="image_25.png" /><figcaption>(a) Chute-du-Diable</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(b) Chute-Savane</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(c) Lac-St-Jean</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(d) Isle-Maligne</figcaption></figure>

Figure 5.8 Water discharges and reservoir levels for the case June 2011, 5 stages, 16 scenarios.

Figure 5.8 pictures June 2011 data set with 5 stages and 16 scenarios. Solutions obtained from the scenario tree method and the median scenario are quite similar. Also note that when a method turbines more water, it is penalized accordingly so it is not advantaged. The absolute difference between the volume, at the end of the 31 day horizon, between the solution obtained with the stochastic model and the solution obtained with the deterministic model is calculated. This volume is then transformed into energy, then added to the method which has a higher end volume, since it is disadvantaged, given the other method processed more water throughout the 31 day planning horizon.

Figure 5.9 also illustrates the June 2011 data set with 7 stages and 16 scenarios. Again, results are very similar.

Without any surprise, the numerical experiments reveal that the solutions to the cases with more stages are closer to the operational ones because the hydropower system operation is
<figure><img src="image_29.png" /><figcaption>(a) Chute-du-Diable</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(b) Chute-Savane</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>(c) Lac-St-Jean</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>(d) Isle-Maligne</figcaption></figure>

<caption>Figure 5.9 Water discharges and reservoir levels for the case June 2011, 7 stages, 16 scenarios.</caption>

more realistic. For example, Figures 5.8 and 5.9 show that the solutions with 5 and 7 stages lead to a similar improvement, but the implementation with 7 stages is preferable. Figures 5.9a, 5.9b and 5.9d present reservoir volumes that are more stable than Figures 5.8a, 5.8b and 5.8d.

The October data set is the only one for which the median scenario produces more energy for all number of stages. The interest of a stochastic method is to account for uncertainty in the future. As we compare our method with the median scenario, if the actual realization of the inflows is close to the median scenario, the stochastic solution will not produce more energy, as the median scenario depicts correctly the future. In practice, this may happen during the fall period, for example when low variability exists in the weather and storms have less chances of developing. This can be seen on Figure 5.10. Each subfigure corresponds to a reservoir. Figure 5.10a is Chute-du-Diable. The top figure is the day 1 October forecast and the bottom figure is the day 1 September forecast. For the first 15 days, the October forecast median scenario is very close to the inflow realization and therefore, as we keep the day 1 decision only, the median scenario produces more energy. The other subfigures