A Modelica-based framework for
modeling and optimization of
microgrids
Ali Akhlaghi
Master of Science Thesis
KTH School of Industrial Engineering and Management Energy Technology TRITA-ITM-EX 2019:543
Division of Heat and Power Technology SE-100 44 STOCKHOLM
A Modelica-based framework for
mod-eling and optimization of microgrids
Ali Akhlaghi
Approved
Examiner
Supervisors
Date: 26/8/2019
Anders Malmquist
Stéphane Velut
Oskar Nilsson
Håkan Runvik
Anders Malmquist
Abstract
Microgrids have lately drawn much attention due to their considerable financial benefits and the increasing concerns about environmental issues. A solution that can address dif-ferent engineering problems - from design to operation - is desired for practical reasons and to ensure consistency of the analyses. In this thesis, the capabilities of a Modelica-based framework is investigated for various microgrid optimization problems. Various sizing and scheduling problems are successfully formulated and optimized using nonlinear and physical component models, covering both electrical and thermal domains. Another focus of the thesis is to test the optimization platform when varying the problem formu-lation; performance and robustness tests have been performed with different boundary conditions and system setups. The results show that the technology can effectively han-dle complex scheduling strategies such as Model Predictive Control and Demand Charge Management. In sizing problems, although the platform can efficiently size the compo-nents while simultaneously solving for the economical load dispatch for short horizons (weekly or monthly), the implemented approach would require adaptations to become efficient on longer horizons (yearly).
Contents
List of Figures 6 List of Tables 10 List of Symbols 11 1 Introduction 14 1.1 State of art . . . 16 1.2 Objectives . . . 17 1.3 Methodology . . . 18 1.4 Structure of thesis . . . 19 2 Background 20 2.1 Hierarchical acausal modeling in Modelica . . . 202.1.1 Software and tools . . . 24
2.2 Differential algebraic equations . . . 25
2.3 Dynamic optimization . . . 26
2.3.1 Optimica . . . 26
2.3.2 Numerical methods . . . 27
2.5 Summary . . . 29
3 Dynamic modeling of microgrid 31 3.1 Electrical system . . . 31
3.1.1 Irradiance and weather data . . . 32
3.1.2 Incident irradiance on PV arrays . . . 34
3.1.3 PV module DC output . . . 36 3.1.4 Battery . . . 38 3.1.5 AC/DC converters . . . 40 3.1.6 Diesel generator . . . 41 3.1.7 Electrical Network . . . 41 3.2 Thermal system . . . 42 3.2.1 Thermal tank . . . 42 3.2.2 Heat exchanger . . . 43
3.2.3 Pumps and volumes . . . 44
3.2.4 Cogeneration plant . . . 44
3.3 Modeling in Dymola . . . 45
3.4 Model input data . . . 47
4 Control methods 51 4.1 Linear quadratic regulator . . . 51
4.2 Model predictive control . . . 52
5 Optimization problem suite 56 5.1 Sizing problem . . . 57
CONTENTS
5.1.1 Sizing the number of PV modules . . . 58
5.1.2 Sizing of the battery capacity . . . 59
5.1.3 Sizing of the diesel generator . . . 60
5.2 Operation scheduling . . . 62
5.2.1 Electrical dispatch problem with MPC and DCM . . . 62
5.2.2 Scheduling of thermal system . . . 64
5.2.3 Scheduling of thermal and electrical combination . . . 65
6 Results 67 6.1 Sizing PV modules . . . 67
6.2 Sizing battery capacity . . . 69
6.3 Sizing diesel generator . . . 71
6.4 Scheduling electrical system . . . 73
6.5 Scheduling thermal system . . . 76
6.6 Scheduling combined electrical and thermal system . . . 78
7 Discussions 81 7.1 Robustness . . . 82
7.2 Performance . . . 83
7.3 Benchmarking of the framework . . . 85
7.4 Sustainability analysis . . . 88
8 Conclusion and future work 90
List of Figures
2.1 Connection diagram of a simple circuit model [20] . . . 21
2.2 Causal block-oriented modeling of the simple circuit [20] . . . 21
2.3 Modeling of simple circuit in Dymola environment . . . 24
2.4 Optimization process of the thesis . . . 30
3.1 Schematic of grid-connected solar PV with battery . . . 32
3.2 The angular position of the sun relative to earth . . . 33
3.3 Sun zenith, elevation and azimuth angle . . . 34
3.4 Equivalent circuit of solar cell . . . 36
3.5 Model of battery efficiency respect to battery power flow . . . 40
3.6 PV production, thermal and electrical load profiles . . . 49
3.7 Electrical and thermal tariffs . . . 49
3.8 The efficiency curves of cogeneration plant based on percentage of load . . 50
4.1 Principle of model predictive control . . . 53
4.2 MPC diagram in OCT . . . 55
5.1 Electrical microgrid with battery and PV in Dymola . . . 58
LIST OF FIGURES 5.3 Modeling of thermal system connected to the network and cogeneration
plant in Dymola . . . 65
5.4 Modeling of combined electrical and thermal system in Dymola . . . 66
6.1 Verification of optimization assumptions . . . 68
6.2 Result of PV sizing problem compare to ELD . . . 68
6.3 Total cost reduction in SPV compare to ELD and Base scenarios. The optimization reduces the PV panel size and therefore the capital cost. The increase in grid power import is however smaller than the savings on the PV cost. . . 69
6.4 Result of battery capacity sizing problem compare to ELD . . . 70
6.5 Total cost reduction in SBC compare to ELD. The battery sizing opti-mization leads in this specific problem to a minimal reduction in the total cost. This reduction is the result of a higher battery capacity and a better management of the grid with less power purchase. . . 70
6.6 The influence of the fuel consumption limit (5000 l) on different scenarios. The plateaus in the profiles are related to stops in the engine. When the fuel budget constraint is activated (SGFB and FB), the maximum fuel amount is spent at the end of the optimization. . . 71
6.7 The sizing optimization (SGFB) leads to a smaller engine that runs at higher loads and more frequently. In the FB scenario (nominal engine size with fuel budget), the engine is used sporadically compared to the unconstrained case (ELDG). . . 72
6.8 The efficiency of the optimally sized engine (SGFB) is higher compared to the other cases (FB and ELDG). This is a direct consequence of the smaller size that the optimization results in. . . 72
6.9 Total cost reduction in SGFB compare to FB and ELDG scenarios . . . 73
6.10 The effect of time horizon of decision variable battery power. Smaller time horizons means shorter forecasts and more aggressive changes in the battery powers when new measurements become available. MPC with 12h horizon is almost identical to the offline optimization over 96h. . . 74
6.11 The effect of DCM on power from grid. The peak power consumed from the grid (around 12h) is reduced by adding the demand charge management strategy in MPC control . . . 75 6.12 The effect of DCM on battery power dispatch. The power from the grid
is reduced by 20% around 12h when formulating the optimization as a demand charge reduction problem (DCM). . . 75 6.13 Cogeneration power dispatch in thermal system. The cogeneration load is
cyclic and follows the thermal grid tariffs. The heat supplied by the grid and the cogeneration plant exceeds the demand during low grid prices (ex: 20-30h), the excess being stored in the tank. . . 76 6.14 Result of power dispatch in storage tank. The link between tariff and use
of storage is clear. The tank is charged at low heat price periods (0-8h, 20-30h) and discahrged at high prices. . . 77 6.15 The thermal storage system flow trajectory . . . 77 6.16 Customer inlet flow trajectory. The temperature of the flow does not
vi-olate the maximum and minimum constraints. The controller keeps the temperature close to lower bound during high heat price periods (8-21h, 32-45h) by increasing the mass flow rate . . . 78 6.17 Thermal production of cogeneration plant. In the combined system (TELD),
cogeneration plant works full-load whole period in contrast to thermal sys-tem (TLD) in which the plant only operates during high heat prices (8-21h, 32-45h). . . 79 6.18 Thermal power flow. While the usage of network decreases in both scenarios
during high heat price periods (8-21h, 32-45h) (Figure 6.18a), the maximum capacity of thermal storage is used to fulfill the thermal demand (Figure 6.18b). . . 79 6.19 SoC of battery profile for TELD and ELD . . . 80 6.20 Operational cost of ELD, TLD and TELD scenarios. . . 80
LIST OF FIGURES 7.1 SoC trajectory with different efficiency. eff = 0.9 represents round trip
efficiency of 0.81. The battery model can operate efficiently based on the
different boundary values for its efficiency. . . 82
7.2 The effect of price on the PV size modifier. Increasing of the price decreases the optimal number of PV panels. . . 83
7.3 SoC of the battery for both frameworks . . . 86
7.4 Battery power for 6 consecutive days . . . 87
List of Tables
3.1 Port definition . . . 45
3.2 Electrical components summery . . . 46
3.3 Thermal components summery . . . 47
3.4 Parameters inputs . . . 48
7.1 Performance of different scenarios . . . 84
7.2 Effect of number of nodes on performance . . . 85
7.3 Comparison of the MPC controllers . . . 87
List of Symbols
Abbreviations
CEC California Energy Commission DAE Differential Algebraic Equation DCM Demand Charge Management DHI Diffused Horizontal Irradiance DNI Direct Normal Irradiance DOP Dynamic Optimization Problem DP Dynamic Programming
ELD Electrical Load Dispatch
ELDG Electrical Load Dispatch with Generator EoT Error of Time
FB Fuel Budget
FMI Functional Mock-up Interface FMU Functional Mock-up Unit
GAMS General Algebraic Modeling System GHI Global Horizontal irradiance
IEA International Energy Agency IP Integer Programming
IPOPT Interior Point OPTimization KKT Karush-Kuhn-Tucker
LMTD Logarithmic Mean Temperature Difference LP Linear Programming
LQR Linear Quadratic Control
LST Local Solar Time LVH Lower Heating Value MGDT Micro Grid Design Tool
MINLP Mixed Integer Nonlinear Programming NOCT Nominal Operating Cell Temperature NTU Number of Transfer Units
OCT Optimica Compiler Toolkit ODE Ordinary Differential Equations PV Photovoltaic
SBC Sizing Battery Capacity
SGFB Sizing Generator with Fuel Budget SoC State of Charge
SPV Sizing PV
SRC Standard Radiation Condition STC Standard Test Condition TC Time Correction
TELD Thermal and Electrical Load Dispatch TLD Thermal Load Dispatch
UTC Universal Time Coordinated Greek letters
α Temperature coefficient β Tilt angle
cos ϑ Share of Gbeonpanels
δ Declination angle Effectiveness
η Efficiency γ Elevation angle ω Cost coefficient ψ Azimuth angle ρ Density τ Transmittance θ Zenith angle υ Algebraic variables ε Tariff price ϕ Latitude ξ Differential variables Latin letters A Surface area
a Absorptance of the PV array B Angle correspond to day of the year C heat capacity of mass flow rate
c Speed of light in a vacuum inertial frame cap Capacity
Cp Specific Heat Capacity d Day of the year
F Diesel consumption coefficient f Derating factor
G Irradiance on the PV array h Enthalpy
HRA Hour angle I Current
Io PV cell diode reverse saturation current K Coefficient of heat loss
L Load l Longitude
m Mass flow N Number
ni Modified nonideality factor op Operational cost
P Electrical power p Parameters Q Heat power
R Regulating decision variables T Temperature
t Time
U Coefficient of heat transfer u Decision variables V Voltage V Volume x State variables Y PV rated capacity Subscripts amb Ambient b Battery be Beam C Cell c Charging c Customer
cog Cogeneration plant d Discharging di Diffuse el Electrical gen Diesel generator h Horizon
he Heat Exchanger in Inlet
LIST OF SYMBOLS
inv Inverter L Light m PV module max Maximum Value min Minimum Value mp Maximum power neg Negative
N et Electrical or thermal network out Outlet p PV panels pos Positive pp Peak power s Sun se Series sh Shaving su Supposed t Tank th Thermal tot Total u Update
Chapter 1
Introduction
With the advent of cost-effective renewable technologies and due to problems related to climate change, reforming the energy industry has become an international priority [1]. According to the International Energy Agency (IEA), renewable energy production will reach 31% of the world’s total power generation by 2035 [2]. The transition to decentralized electricity generation is a crucial step for economic and environmentally friendly energy production. Traditionally, electricity is generated at large-scale power plants and distributed to the consumers via a network of high and low-voltage transmission and distribution lines. Shifting to a distributed generation, often intermittent such as small-scale photovoltaic (PV) and micro-turbines, usually located closer to the consumers, has a number of potential benefits, including reduced peak power consumption, increased power quality and reduced transmission loss. In fact, energy production is becoming more decentralized and energy consumers are becoming so called prosumers that also feed energy back into the grid and store energy locally.
The important challenge of decentralized renewable generation is dealing with the increase in non-dispatchable energy sources that makes it more difficult to balance energy supply and demand. Dispatchable energy sources refer to those sources with a load that can be controlled and predicted, for instance fuel-based energy systems like diesel generator, micro gas turbine and fuel cell. In contrast, non-dispatchable energy sources, which include intermittent energy sources such as wind, solar, and wave energy, impose more volatility to the power grid due to their unforeseen nature. The large-scale deployment of distributed solar PV panels is, as an example, very challenging because of the dependency of the solar energy with geographical area and sun diurnal motion. PV output generally has a peak production during the midday sun and then subsides in the late afternoon
as the sun sets. Without storage, the electricity generated by the midday sun cannot be used to meet the residential electricity demand in the late afternoon and evening. According to [3], the self-consumption rate of solar panels -the ratio between the load and the available solar energy PVs- was in 2014 between 17% and 50%. This large variation in solar power and the mismatch in times of production and demand show that, on one hand, weather and load demand forecasting and, on the other hand, thermal and electrical storage are crucial components to improve energy management systems and to increase the self-consumption rate of microgrids.
Typically, a microgrid comprises four units including, storage units, flexible energy loads, dispatchable energy generators such as distributed generators and non-dispatchable energy generators for instance photovoltaics and wind turbines. One of the main motivations behind using microgrids is their capability of being controllable in a more decentralized way, which decreases the dependency on centralized coordination and management; as a result, optimization of its operation is crucial. Moreover, due to the considerable capital cost of energy systems especially electricity storage, studies have paid great attention on the optimal sizing of energy components [1]. Hence, introducing a unique framework which can tackle both optimal sizing and dispatch problems would be beneficial.
This project is done as a master thesis in collaboration with the two companies, Schneider Electric and Modelon [4, 5]. Schneider Electric leads the digital transformation of energy management and automation in housing, real estate, data centers, infrastructure and industries. With presence in over 100 countries, it is one of the largest companies in energy management for medium voltage, low voltage uninterrupted power and automation system. It offers integrated solutions that effectively combine energy, automation and software. In this master thesis, Schneider electric in collaboration with the simulation company Modelon is planning to use the language Modelica in order to model and optimize different microgrid problems.
Modelon offers a comprehensive suite of model libraries, built on the Modelica standard, that delivers state-of-the-art system models for a wide range of industries including au-tomotive, aerospace, industrial equipment, and energy and process. Modelon libraries are developed and updated in close cooperation with leading industry partners to reflect the evolving industry needs. Modelon is also developing a comprehensive simulation and optimization web-based tool to provide a user-friendly solutions for diverse industries.
1.1
State of art
Four common problems are identified for an economic feasibility approach to microgrid, including, power mix selection, sizing, siting and scheduling [6]. The first three problems are considered as strategic level problems while the last one focuses on tactical planning level. Power mix selection is the problem of allocating the available resources of the project location to setup the best power system fulfilling the requirements. Sizing deals with the optimal dimensioning of the available equipment. Sitting aims at keeping quality constraints in the power lines. Finally, scheduling focuses on dispatching the production between the available units and mostly aimed at minimizing the operational cost while ful-filling operational and environmental constraints. This thesis tries to show the capability of the Modelica-based framework to address the sizing and scheduling problems.
Introducing a unified tool-chain capable of handling all the relevant problems in both design and operation is a necessity for optimal planning of microgrids. After the selection of power sources, the computation of the optimal components’ size is an important factor that influence both the capital and operational costs of the microgrid, which is the reason for solving the sizing problem in connection with the economic load dispatch problem [6]. Different methods are proposed to solve sizing problem along with economical load dispatch. In [7], the problem is formulated as a mix integer programming (MIP) to minimize the operation cost. Heuristics are also widely used in sizing problem. In [8], as an example, some heuristic optimization techniques are implemented for hybrid renewable energy system sizing such as: PSO, SA, GA and Ant Colony.
Furthermore, in the literature different algorithm and control strategies are proposed to tackle the scheduling problem. Fuel consumption rate is minimized while obeying to fulfill the local energy demand and provide extra backup power by implementing MINLP in [9]. Studies [9, 10] , shows depend on problem definition IP and LP are suitable approach. Dynamic programming is used to tackle both types of problem, in [11] for sizing and in [12] for scheduling problem. Since inequality constraints and nonlinearities can be included in microgrid problem definitions, it is an attractive method for optimization formulation of real cases. Moreover, in scheduling problems, optimization algorithms can make large improvements in the operation and also substantially reduce the costs, especially if the demand charge can be accounted for in the optimization. The demand charge is the price for industrial building with high load profiles that has to be paid for the maximum power consumption in off, mid and on-peak periods. The demand for online control systems that can deal with the abrupt changes in renewable energy resources is also an argument for the application of optimal control.
1.2. Objectives Platforms for microgrid management require certain features in order to tackle variety of problems. Currently, applications are widely used to simulate and optimize microgrids such as HOMER [13] and PSCAD [14]. HOMER (Hybrid Optimization of Multiple Energy Resources) is a commercial tool originally developed by the National Renewable Energy Laboratory (NREL). It simulates the system to create a list of feasible configuration which can be sorted by different system-defined criteria to locate the optimal solution. This makes Homer a useful and efficient tool for optimization of the system configuration, component sizing and offline grid operation control. PSCAD is a time-domain based power system simulation tool, which can be applied for microgrid configuration. However, these tools are quite specialized, lacking the needed flexibility of a complete solution, such as changing fidelity levels, customizing optimization formulations, adding highly customized component models and implementing the combination of diverse physical domains [15] while keeping its user-friendly characteristics. For instance, the work required to setup such systems with the generic and powerful modelling language GAMS (General Algebraic Modeling System) [16] is too large and makes such tools impractical for the common user. Modelica is primarily a modeling language that targets the modeling of complex physical systems. A Modelica-based solution is a potential candidate to fill the lacking features needed in microgrid applications, by offering a highly flexible and user-friendly framework. In [15], the authors describe the new venture of a Modelica and web-based solution for microgrid analysis and design. This study is further developed in this thesis by introducing more real cases problems, evaluating the performance in different cases and benchmarking with other tools.
1.2
Objectives
The goal of the project is to investigate the capacity of a Modelica framework in han-dling diverse kinds of microgrid related problems. Premiere focus is to further develop the Modelica-based microgrid library built by Modelon in order to make a robust and comprehensive library being able to handle various optimization problems. Thereafter, the goal is to study different optimization problems that can show the strength and the eventual gaps of the platform. Two types of problems, i.e. sizing and scheduling are cho-sen to assess the framework. Different configuration sets for each problem are considered and solved. The assessment of the framework is completed by analyzing the results from the different optimization problems.
models such as thermal domain in the problem formulation. Most of the current platforms require significant modifications in order to include nonlinear thermal components along with linear and nonlinear electrical ones. However, the Modelica-based tool proposes a flexible platform in which differentiable model regardless of being linear or nonlinear can be easily added, by drag-and-drop, to any existing configuration. It facilitates and accelerates the design process of complex systems including the combination of physical domains.
Furthermore, the Modelica-based tool is benchmarked with other solutions, particularly Schneider Electric micro grid design tool (MGDT) for a more complete assessment of the the platform performance in comparison with current available solutions. The benchmark provides suitable information about the platform and can illustrate its strengths and weaknesses. Computation time, as an example, is a crucial factor for any alternative application that can show the growth capacity of the platform among users.
To put in a nutshell, goals of the thesis can be summarized as follows: • Further development of the microgrid library by Modelon.
• Assessment of the framework through the optimization of different single or multi-domain (electrical and thermal) suites in sizing and scheduling problems.
• Benchmarking of the tool with MGDT.
1.3
Methodology
The focus of the thesis is the optimization of microgrids in the Modelica framework by Modelon and introducing new models for missing components. New components are modelled based on the review of previous studies on modeling of microgrid and common tools such as Homer energy. The implementation of previously tested physical models from other projects or studies can provide the opportunity to verify the results generated by the new approach. However, since the Modelica-based optimization platform demands all models to be continuously differentiable, some modifications of the original models are required to make them compatible with the tool.
Acausal and physical block-oriented modeling is used here to model dynamic systems of connected components in contrast to block-oriented, causal modeling which is the classical algorithm-based approach. The causal approach makes it possible to build hierarchical
1.4. Structure of thesis structure of input-output blocks but it hides the physical reality under the computational structure. In the acausal Modelica tools, the components are coupled using physical connections, which gives the model a structure that is similar to the real system and does not require the user to rewrite the model as sequential and cause-and-effect computational flows. Acausal block-oriented modeling generates differential algebraic equations (DAEs), which are the combination of differential equations and algebraic equations.
The optimization problems are defined using Optimica, an extension of the Modelica language that allows the user to define cost function and constraints. The numerical approach, is a so-called collocation method that discretizes the continuous-time description of the dynamic optimization problem to get a sparse discrete-time nonlinear program. In contrast to typical micro-grid optimization tools like Homer that implement standard control strategies and find the best system configuration by ranking results from massive simulations, the Modelica-based solution finds the optimum using an explicit and gradient based method on all system equations. The approach has been successfully applied to other energy applications like start-up of power plants or production planning of district heating systems [17, 18].
1.4
Structure of thesis
Chapter 2 covers the tools and the methods used in the thesis for the modeling and optimization of the microgrids: the Modelica language, system modelling in Modelica, the optimization method and the optimization process. Chapter 3 is dedicated to the review of the microgrid components used in the optimization problems and describes the details of the models both for the electrical and thermal systems. Chapter 4 is focusing on the control strategies used in the thesis. Chapter 5 includes the configurations of the system to be optimized and the problem formulations. Chapter 6 provides the results from the optimization problems together with their analysis. Chapter 7 makes a discussion on reliability of models and system configuration. Moreover, it contains some performance analysis for different problems and some benchmark with other tools. Finally, Chapter 8 is the conclusion and suggestion for future developments.
Chapter 2
Background
In this chapter, the background material is reviewed with the purpose of introducing the main concepts of nonlinear optimization methods to facilitate the comprehension of the thesis. The thesis consists basically of two main parts: modeling and optimization. The chapter starts with theories and programs for the modeling part in Modelica, then it continues with some theoretical background needed for understanding the optimization process in Optimica Compiler Toolkit (OCT).
2.1
Hierarchical acausal modeling in Modelica
Block-oriented, causal modeling is the classical approach to model dynamic systems. In this approach, the model is written as cause-and-effect flows between input and output signals, which can make it difficult to represent complex systems due to the restructur-ing of all physical equations [19]. Although the correspondrestructur-ing mathematical formalism of causal modeling can be computationally efficient since it results in explicit ordinary differential equations (ODEs), it may be inconvenient for the modeling of systems that are intrinsically described by differential algebraic equations (DAEs). Furthermore, causal modeling prevents model reuse because every model is built for a specific cause-and-effect flow, which is an important drawback for making a comprehensive library of components to be used in various ways.
Figure 2.1, an example of a simple circuit, is presented to illustrate the idea of causal modeling disadvantages (from [20]) and is modeled based on causal block-oriented model-ing depicted as a diagram in Figure 2.2. The model comprises a sinusoidal AC source AC,
2.1. Hierarchical acausal modeling in Modelica two resistors R1 and R2, a capacitor C, an inductor L and ground G. Since the model is causal (Figure 2.2), the signal flow is straightforward and the output can be calculated easily by defining the input. However, it is evident that, comparing the structure of the original physical circuit and the block diagram, the physical topology is lost. Moreover, the resistor blocks are context dependent meaning that, for example, the resistors R1 and R2 have different declarations, which makes the reuse of model library components difficult. Another negative point of this modeling approach is that system models need to be considerably re-built when minor changes or additions are done to the nominal system. Finally, even for this simple example, conversion of physical model to a causal block-oriented model requires a nontrivial analysis. To address those disadvantages, the block-oriented acausal modeling languages such as Modelica can be used.
Figure 2.1: Connection diagram of a simple circuit model [20]
Modelica is primarily a modeling language with the target of modeling heterogeneous physical systems. It provides a framework to describe mathematical models of complex nature or user-defined systems and it can, for instance, be used for dynamic simulation or optimization of multi-domain systems. Furthermore, Modelica is an object-oriented equation-based programming language that simplifies the modeling of complex systems consisting of several sub-components and a hierarchical structure. The four noticeable features of Modelica are [20]:
• Modelica is equation-based in contrary to assignment statements that give a value to a variable. This characteristic ameliorates reuse of classes as there is no need to specify a certain data flow direction; as the result, the models can adapt to several data boundary conditions.
• Modelica is capable to model multi-domains, meaning that a model can be a con-nection of sub-components corresponding to several physical domains such as, e.g., mechanical, electrical, thermodynamic, hydraulic, and control applications.
• Modelica is an acausal object-oriented language with a standardized class format, facilitating reuse and exchange of models.
• Modelica is supported by several software with graphical user interfaces to construct and connect different components. Thus, it is suited as an architectural description language for complex physical systems.
All above mentioned properties make Modelica a powerful and user-friendly technology for modeling complex systems. To illustrate the hierarchical acausal object-oriented modeling in Modelica, the circuit in Figure 2.1 again is modeled, this time as a Modelica class called SimpleCircuit in Listing 2.1. The model is acausal since no signal flow can be detected. The physical structure of system is maintained in the model because by drawing a simple analogy each object in the circuit can be connected to a declared instance in the SimpleCircuit. Connections between objects are specified using connect function in an equation section. Therefore, it provides a framework to reuse the models in the case of changes in the physical structure.
2.1. Hierarchical acausal modeling in Modelica
Listing 2.1: Acausal block-oriented modeling in Modelica
m o d e l S i m p l e C i r c u i t R e s i s t o r R1(R=10); C a p a c i t o r C(C=0 .01); R e s i s t o r R2(R=100); I n d u c t o r L(L=0 .1); V s o u r c e A C AC ; G r o u n d G ; e q u a t i o n c o n n e c t(AC.p , R 1 . p); // C a p a c i t o r c i r c u i t c o n n e c t(R1.n , C.p); c o n n e c t(C.n , A C . n); c o n n e c t(R1.p , R 2 . p); // I n d u c t o r c i r c u i t c o n n e c t(R2.n , L.p); c o n n e c t(L.n , C.n); c o n n e c t(AC.n , G.p); // G r o u n d
The hierarchical structure of modeling in Modelica can be demonstrated in the simple cir-cuit example. The components Resistor, Capacitor, Inductor, VsourceAC and Ground have already been defined in different classes meaning that they may include some sub-components inside and imported to top level for final configuration. For instance, the component Resistor is made from a connection pin named TwoPin and an equation based on ohm’s law which relates the boundaries value, voltage and current, to the intrin-sic parameter of component, ohmic resistor.
Listing 2.2: Resistor model
c l a s s R e s i s t o r " I d e a l e l e c t r i c a l r e s i s t o r " e x t e n d s T w o P i n ; p a r a m e t e r R e a l R(u ni t=" Ohm ") " R e s i s t a n c e "; e q u a t i o n R*i = v ; end R e s i s t o r ;
TwoPin class being extended to the higher-level works as a connection interface for elec-trical components and is defined as follows:
Listing 2.3: Electrical interface model
p a r t i a l c l a s s T w o P i n " S u p e r c l a s s of e l e m e n t s wi t h two e l e c t r i c a l p i n s "
Pin p , n ; V o l t a g e v ; C u r r e n t i ;
v = p.v - n.v ; 0 = p.i + n.i ; i = p.i ;
end T w o P i n ;
Each Modelica class comprises two main parts: the variable declaration section and the equation section. The equation section allows the user to model a component by hybrid differential-algebraic equations. This means that differential equations, algebraic con-straints between variables, as well as discrete-time events can be combined and simulated using Modelica. The Modelica components are connected with well-defined connectors to build a more complex component. Well-defined connectors and interfaces allow com-patible components to be connected to each other, and makes it easy to build complex architectures.
2.1.1
Software and tools
The software used in this thesis to build Modelica models is Dymola [21] which is a commercial program offering both a modeling environment and a Modelica translator. Models can be created both in the text view, by editing the Modelica code, or in the diagram view, by connecting components from the graphical interface. Figure 2.3 shows the modeling environment of Dymola for the simple circuit model(Figure 2.1). It should be noted that several component libraries, either open-source or commercial, are available for Dymola. In this project, for instance, the commercial microgrid library designed by Modelon is used for modeling the physcial systems.
2.2. Differential algebraic equations When modeling is completed, to exchange of components between users and tools, the component can be packaged in a Functional Mock-up Unit (FMU), a binary file complying with Functional Mock-up Interface (FMI) standards [22]. A tool that implements the FMI standard can use the model, alone or as a part of a system, to solve a simulation problem. The FMU can also include a solver, which is then referred as FMU for Co-Simulation (FMU-CS). In the thesis, Optimica compiler toolkit (OCT), Modelon’s Modelica and FMI computational platform, is used to compile and load FMUs[23]. Moreover, it provides an environment to solve dynamic and steady-state simulation and optimization problems. From the next section, the theoretical background for optimization in OCT is covered.
2.2
Differential algebraic equations
After restructuring physical models to create a clear signal flow between inputs and out-puts in causal, block-oriented modeling approach, generally the mathematical formalism is explicit ordinary differential equations (ODE) of the form of
˙x(t) = f (t, x(t), u(t), p) (2.1)
Where t is time, x is the state variable, u is decision input and p represents system parameters. Most ODEs are solvable numerically with common ODE solvers. On the other hand, acausal, hierarchical modeling is a generalization of implicit ODEs in which the differential variables are defined implicitly in the format of differential algebraic equations DAEs. Therefore, the equations set consists of both differential equations and algebraic equations. Implicit DAEs are of the form of
Φ(t, ˙ξ(t), ξ(t), υ(t), u(t), p) = 0 (2.2) Where ξ is the differential variable and υ is the algebraic variable. Based on the thesis focus on microgrids, in both thermal and electrical domains differential variables are explicitly formulated which leads to another form of DAEs called the semi-explicit form
˙
ξ(t) = f (t, ξ(t), υ(t), u(t), p) (2.3a)
g(t, ξ(t), υ(t), u(t), p) = 0 (2.3b)
2.3
Dynamic optimization
Dynamic optimization problems (DOP) are generally formulated in the form of combina-tion of semi-explicit DAEs along with constraints and objective funccombina-tion [24].
minimize φ(tf, ξ(t)) + Z tf t0 L(ξ(t), υ(t), u(t))dt, (2.4a) with respect to ξ : [t0, tf] −→ Rnξ, υ : [t0, tf] −→ Rnυ, u : [t0, tf] −→ Rnu, tf ∈ R, p ∈ Rnp, subject to ξ(t) = f (t, ξ(t), υ(t), u(t), p),˙ (2.4b) g(t, ξ(t), υ(t), u(t), p) = 0, (2.4c) ξ(t0) = ξ0, (2.4d) h(t, ˙ξ(t), ξ(t), υ(t), u(t), p) ≤ 0, (2.4e) H(t, ξ(t), υ(t), u(t), p) ≤ 0, (2.4f) ∀t ∈ [t0, tf]
where 2.4a is the objective consisting the Mayer term (final penalty) φ and Lagrange integrand L, 2.4b and 2.4c are the explicit DAE, 2.4d refers state initial condition, 2.4e is inequality constraint, and 2.4f is the terminal constraint. u, tf, and p are the degrees
of freedom determining ξ and υ by the system dynamics 2.4b and 2.4c.
Since Modelica is mainly designed for simulation-based analysis, it is unable to support optimization formulations; as a result, Modelica is only able to encode 2.4b to 2.4d which are DAEs and initial conditions. To address the DOP formulation, Optimica, an extension of the Modelica language, can be used.
2.3.1
Optimica
Optimica is an extension of the Modelica language used for dynamic optimization of Modelica models [23]. It contains a specialized class for optimization in which an opti-mization problem can be defined. The class contains the four new attributes objective, ObjectiveIntegrand, startTime and finalTime. The first and second encode φ, L in 2.4a and the two remaining encode time horizon. Furthermore, a new section called constraint is included to describe the 2.4e and 2.4f sections, which is similar to equa-tion secequa-tion in Modelica (Listing 2.2) with the ability to contain inequalities and equality constrains. The constraints can be point-wise or valid over the optimization period.
2.3. Dynamic optimization In Modelica, each variable can include two attributes called min and max which are the lower and upper bounds. These characteristics are only for representing the range of each variable and is not active to change the result during the simulation. However, by extending the model to Optimica for creating an optimization problem, these attributes become active and work as additional constraints in the DOP
In Optimica, optimization of one or more parameters along with decision variables sizing can easily be defined without requiring changes of the model. For instance, sizing problem based on economic load dispatch is to optimize DOP formalism 2.4 with two different kinds of optimization variables including, one or more decision variables u and one or more parameters in the system p. On one hand, decision variables can be easily inserted to the Dymola model by defining them as input. During simulation, also they are controlled with a feedback controller. On the other hand, in order to reuse the model, parameters must be defined inside the Dymola model with their corresponding values. Then, in Optimica, each parameter can be converted to a degree of freedom by activating the free attribute of the parameter. Therefore, for shifting from a scheduling problem to a design problem in microgrids, there is no need to change the model and it can be easily done during optimization formulation.
2.3.2
Numerical methods
In this thesis, the solution of the DOP is (approximately) found with the direct collocation numerical method which is the way of OCT for handling dynamic optimization of DAEs. Since the focus of the thesis is more on microgrid systems, a brief description of the method is adequate for understanding of the following sections.
Direct methods in contrast to indirect methods, first discretize the system dynamics which reduces the DOP to a nonlinear program, and then impose the optimal condition. The discretization is done by the collocation method in which the dynamic model variable pro-files are approximated by piece-wise polynomials. The elements divide the optimization time horizon into smaller sections and the state profiles contain several Radau collocation points inside an element. The polynomials are chosen to match the values and derivatives of the states at each collocation point[25]. For implementation of the direct collocation method, OCT uses the automatic differentiation tool CasADi which is a symbolic frame-work for nonlinear optimization and algorithmic differentiation [26].
2.4
Nonlinear programming
In previous section, direct methods transfer the indimensional DOP into a finite-dimensional nonlinear program (NLP) with the general form of
minimize f (x), (2.5a)
with respect to x ∈ Rnx,
subject to xL≤ x ≤ xU, (2.5b)
g(x) = 0, (2.5c)
h(x) ≤ 0, (2.5d)
The NLP solvers require that the model equations are twice continuously differentiable with respect to all variables [27]. The most conventional methods for solving NLPs are either active-set sequential quadratic programming or interior-point methods which are approximating the NLP by an equality-constrained NLP. Since OCT uses the IPOPT package, which refers to Interior Point OPTimization, in the following its basics are ex-plained [28].
2.4.1
Interior-point methods
A conventional step is to transfer the inequality constraints by introducing a slack variables to obtain minimize f (x), (2.6a) with respect to x ∈ Rnx , s ∈ Rns, subject to xL≤ x ≤ xU, (2.6b) g(x) = 0, (2.6c) h(x) − s = 0, (2.6d) s ≤ 0. (2.6e)
For better illustration, a simpler NLP is considered as below
minimize f (x), (2.7a)
with respect to x ∈ Rnx,
subject to g(x) = 0, (2.7b)
2.5. Summary A common variant of interior-point methods uses logarithmic functions to achieve the equality constrained NLP [25], minimize f (x) − µ n X t=1 ln(xi), (2.8a) with respect to x ∈ Rnx, subject to g(x) = 0, (2.8b) (2.8c) where µ is a barrier parameter greater than zero and by tending to zero, the optimization problem is solved. To solve the NLP, it is important to define the Lagrangian function
L(x, λ) := f (x) − µ
nx
X
t=1
ln(xi) + g(x)Tλ. (2.9a)
The optimal solution must satisfy Karush-Kuhn-Tucker (KKT) being given by [25], ∇xL(x, λ) := ∇f (x) − µx−1+ ∇g(x)Tλ, (2.10a)
g(x) = 0. (2.10b)
2.5
Summary
The process of optimization of a system is shown in Figure 2.4. First, the system is modeled in Dymola then OCT is used to compile the Modelica generated model to FMU and loaded. At this point, the framework is ready to simulate the model with customized settings (e.g. solver algorithm and tolerance). The main purpose of this simulation is to compute a suitable set of initial trajectories for each variable. It can considerably improve the performance of optimization, but this is often necessary as our problems are not convex, so it is required to start close to the solution. In the next step, the opti-mization formulation is defined in Optimica and is combined with the Modelica dynamic model to create the DOP. Next, by using the direct collocation method with automatic differentiation tool CasADi, the DOP is transfered to NLP and ready to invoke an opti-mization algorithm. Finally, the optiopti-mization algorithm of OCT for dynamic optiopti-mization (IPOPT) is implemented to solve the resulting NLP.
Chapter 3
Dynamic modeling of microgrid
In the thesis, the final objective of the system design is to include both thermal and electrical components in a microgrid. Since the modeling is based on an acausal block-oriented approach, each static or dynamic component must be modeled in a class which represents its physical behaviour. As a result, first, all components needed for a microgrid must be included in a library afterward, every class can be instantiated based on its microgrid setup.
In this chapter, the physical models that are the base for all components used in opti-mization and simulation problems are described. First, electrical components are pre-sented comprising the calculation of solar power output, battery, weather calculation and inverter. Next, thermal components are described including, heat exchanger, tank, co-generation plant, thermal network, pipe volumes and customer with both thermal and electrical demand.
3.1
Electrical system
A photovoltaic system is a power system in which solar energy is converted to electricity by means of photovoltaic. Figure 3.1 shows the schematic of a typical grid-connected PV system with an AC bus connected battery. It consists of several components, including solar panels to absorb and convert sunlight to electricity, inverters to change the electricity current from DC to AC, rectifier to do the opposite of the inverter, which means to alter the electricity current from AC to DC and a battery to improve the overall performance of the system by storing and managing of excess electricity. In the following of this section,
the electrical components models are described in detail.
Figure 3.1: Schematic of grid-connected solar PV with battery
3.1.1
Irradiance and weather data
Since PV arrays generate electricity from the sunlight, it is important to calculate the arrays’ incident solar radiation at each time step from its weather data. Generally, weather data providers have utilized three parameters, including global horizontal irradiance (GHI) direct normal irradiance (DNI) and diffused horizontal irradiance (DHI) to describe the solar irradiance at each location. GHI is the total amount of irradiance received from above by a surface horizontal to the ground. DNI is the power of radiation coming directly from the sun onto a normal unit surface. DHI is the radiation that comes from all parts of the sky, except the sun. It is worth mentioning that it is measured by a horizontal detector (pyranometer) with a mobile shadow ring that blocks the radiation emitted directly from the sun [29]. Besides the irradiance data, ambient temperature (Tamb) is required for more
3.1. Electrical system
Sun orientation and position
In order to calculate the solar incident, first the position of the sun has to be calculated. In the following, the method of Duffie and Beckman ([30]) is used to address the problem. The position of the sun relative to a plane, can be described in terms of several angles. For calculating them, initially from the Equation (3.1) universal time coordinated (UT C) must be converted into the local solar time (LST ). Hence,
LST = U T C + T C
60 , (3.1)
where T C is time correction factor which is calculated from Equation (3.2).
T C = 4l + EoT. (3.2)
In the Equation (3.2), l is the longitude and error of time (EoT ) which is obtained from Equation (3.3).
EoT = 9.87 sin 2B − 7.53 cos B − 1.5 sin B (3.3) By considering, d as the number of the day in the year (d=1 corresponds to January, 1st and d=81 to spring equinox i.e. March, 22nd ), B is computed from the equation below,
B = 360
365(d − 81). (3.4)
Figure 3.2 shows the declination angle (δ) which is the angular position of the sun at solar noon with respect to the plane of the equator, north positive. It can be found from the approximate equation of Cooper (1969) [31],
δ = 23.45 sin 360
365(d − 81)
(3.5)
Figure 3.3 illustrates the solar elevation angle γs, solar zenith angle θs and solar azimuth
angle ψs. They obey the approximate relations [30]:
γs = arcsin(sin δ sin ϕ + cos δ cos ϕ cos HRA) (3.6a)
θs = 90◦− γs (3.6b)
ψs=
arccos(sin δ cos ϕ−cos δ sin ϕ cos HRAcos γ
s ), if HRA ≤ 0
360 − arccos(sin δ cos ϕ−cos δ sin ϕ cos HRAcos γ
s ), otherwise
(3.6c) where ϕ is the latitude and HRA is the hour angle derived from below:
HRA = 15(LST − 12). (3.7)
Figure 3.3: Sun zenith, elevation and azimuth angle
3.1.2
Incident irradiance on PV arrays
Solar arrays must be oriented in the right direction to achieve as much sunlight as possible during a day. Nowadays, three possible mounting alternatives, including the fixed, one-axis tracking, and two-one-axis tracking are implemented. Selecting the best orientation and position is highly related to the project and its location due to the fact that both cost and efficiency of the system will increase as we go from choosing the fixed to two-axis tracking [32].Furthermore, some limitations may prevent utilizing one or two of options. For instance, it is preferable to have fixed mounting on gable roof houses while two-axis and one-axis are more utilized in solar farms or flat roof houses.
3.1. Electrical system Different mounting options only affect tilt βp and azimuth ψp angle of panels in the model.
They are constant for the fixed mounting, βp is constant and ψp is equal to solar azimuth
angle at each timestep for one-axis tracking and βp and ψp are respectively, zenith and
azimuth angle of sun for two-axis tracking at each timestep. Fixed : βp = β0, ψp = ψ0
One-axis tracking : βp = β0, ψp = ψs
Two-axis tracking : βp = θs, ψp = ψs
(3.8)
The total solar irradiance incident on the PV arrays (Gtot) at each time step is calculated
by sum of diffuse irradiance (Gdi) and beam irradiance (Gbe),
Gtot = Gbe+ Gdi (3.9)
At each time step the beam solar irradiance on the panels surface is linked to the direct normal irradiance (DNI) being available from weather data:
Gbe = DN I cos ϑ (3.10)
where cos ϑ is the scalar product between the unit vector pointing towards the sun and the unit vector normal to the surface,
cos ϑ = sin θscos ψs sin θssin ψs cos θs . sin βpcos ψp sin βpsin ψp cos βp. (3.11)
Because cos ϑ is negative when the sun is behind the surface, equation(3.10) is modified to
Gbe = max(DN I cos ϑ, 0). (3.12)
It is assumed that the diffuse radiation from the sky is isentropic; as a result, the share of diffuse irradiance collected by the surface of panels from the total diffuse horizontal irradiance (DHI) is (the albedo of the ground is neglected)[33]:
Gdi = DHI
1 + cos βp
3.1.3
PV module DC output
Solar PV cell DC output has been calculated in various approaches with different levels of accuracy from the PV panel data sheet. Figure 3.4 is the single diode equivalent circuit of a PV cell where VC, IC, Rse, Rsh, Io and IL represent output current of cell, terminal
cell voltage, series resistance, shunt resistance, diode reverse saturation current and light current, respectively.
Figure 3.4: Equivalent circuit of solar cell
Equation (3.14) shows the electrical nonlinear relation between current and voltage of a cell based on its characteristics. The power output of each cell (PC) is the product of
voltage and current of each cell. The calculation of power output from equation (3.14) is not straightforward because most manufacturer data sheets generally include only a standard set of parameters that define several key points on a module’s I–V curve at stan-dard radiation condition (SRC), as well as other information about the module structure and wiring that are insufficient to calculate power output of a cell. Therefore, the level of accuracy of PV power calculation highly depends on how the effects of parameters in Figure 3.4 are considered.
IC = IL− Io(exp VC+ ICRse niC − 1) − VC + ICRse Rsh (3.14) One of the most accurate methods for calculation of PV DC output is presented by W.De Soto et al [34]. The five-parameter model uses data provided by the manufacturer, absorbed solar radiation and cell temperature together with semi-empirical equations, to predict the current–voltage curve. Moreover, Dobos has described an improved algorithm for calculating the six parameters required by the California Energy Commission (CEC) PV Calculator module model [35]. However, due to the goal of thesis which is the study of high-level system performance, to avoid complexity the simple model of Homer energy for solar cells is selected [36]. This method only calculates the power flow of modules and set the terminal voltage as the nominal value. Hence, the PV module power output (Pm)
3.1. Electrical system is calculated from below,
Pm = Ymfm(
Gtot
Gtot,ST C
)(1 + αm(TC − TC,ST C)) (3.15)
where Ym is the rated capacity of PV module under standard test condition(STC), fm
is module derating factor, Gtot,ST C is the solar irradiance on the PV array under STC,
αm is the temperature coefficient of power, TC is the PV cell temperature in the current
timestep and TC,ST C is cell temperature under STC and total power output of PV arrays
is Pm times number of panels (Np).
Pp = NpPm (3.16)
All the above-mentioned parameters can be obtained either from weather data or from the solar cell data sheet except the PV cell temperature. The PV cell temperature is the temperature of the surface of the PV array which is same as the ambient temperature at night. However, due to the thermalization process in PV solar cells when they have been exposed to the solar radiation, it can exceed the ambient temperature by 30 degrees or more.
By considering an energy balance for the PV array which is the solar energy absorbed by the PV array on one hand and the electrical output plus the heat transfer to the surroundings on the other hand [30],
τsasGtot = ηe,pGtot+ Up(TC − Tamb) (3.17)
where τs is the solar transmittance of any cover over the PV array, as is the solar
absorp-tance of the PV array, ηe,p is the electrical conversion efficiency of the PV array and Up is
the coefficient of heat transfer to the surroundings. The equation 3.17 is solved for cell temperature to yield:
TC = Tamb+ Gtot τsas Up 1 − ηe,p τsas . (3.18)
It is difficult to measure the value of τsas
Up directly, so instead manufacturers report the
nominal operating cell temperature (NOCT), which is defined as the cell temperature that results at an incident radiation of 0.8 kW/m2, an ambient temperature of 20◦C, and no
for τsas
Up , yielding the following equation,
τsas
Up
= TC,N OCT − Tamb,N OCT
Gtot,N OCT (3.19)
where TC,N OCT is the nominal operating cell temperature, Tamb,N OCT is NOCT ambient
temperature(20◦C) and G
tot,N OCT is NOCT irradiance(0.8 kW/m2). By the assumption
of constant τsas
Up [37], the equation 3.19 can be substituted to the equation 3.18 leading to:
TC = Tamb+ Gtot(
TC,N OCT − Tamb,N OCT
Gtot,N OCT
)(1 − ηe,p τsas
) (3.20)
Duffie and Beckman (1991) assumes a value of 0.9 for τsas in the equation above [30].
Because the term ηe,p
τsas is small compared to unity, this assumption does not introduce
a significant error. The second assumption is that the PV array always operates at its maximum power point, as it does when controlled by a maximum power point tracker [37] which means the cell efficiency(ηe,p) is always equal to the maximum power point
efficiency(ηmp,p). By considering linear relation between (ηmp,p) and cell temperature, it
can be written,
ηmp,p= ηmp,p,ST C(1 + αm(TC − TC,ST C)) (3.21)
where ηmp,p,ST C and TC,ST C are the maximum power point efficiency and the cell
temper-ature under STC, respectively. The tempertemper-ature coefficient of power is normally negative, meaning that the efficiency of the PV array decreases with increasing cell temperature. Finally, the cell temperature can be calculated by the set of equations (3.17-3.21). By considering the explicit method used for calculation of PV power output based on location and weather data, adding the calculation of PV system from input to the opti-mization framework not only does not bring any advantages for the optiopti-mization problem, but also due to non-differentiable functions for variable such as day of year, it can deteri-orate the performance of the optimization. Consequently, the simulation of PV output is done before the optimization and the value is used as an external input for the problem.
3.1.4
Battery
Battery implementation in a microgrid can have several purposes such as increasing the self-consumption rate and independence from the grid, reduce peak power consumption and ameliorating power quality [38]. In the thesis, the battery is used to first include dynamics in the electrical part of microgrid and then to check the ability for peak shaving application.
3.1. Electrical system A battery adds dynamics to the electrical model because its state of charge(SoC) depends on previous time step value and the total charging or discharging flow. Hence,
SoC(i) = SoC(i − 1) + Z ti
ti−1
ηbIbdt (3.22)
where t is time, Ib is current to or from battery, ηb is efficiency of battery and i is the
time step. Since the PV is modeled based on power flow, the equation (3.22) is rewritten based on power flow,
d(SoC) dt = ηb
Pb
capb
(3.23) where Pb is battery power flow and capb is battery capacity. In the dynamic modeling of
the battery, efficiency must be implemented differently during charging and discharging due to the direction of flow. When the battery is in charging mode, the actual power affecting the SoC is less than the power going to the battery unit. On the other hand, in discharging mode, the power extracted from SoC is more than real power flow from the battery unit. Hence,
ηb = ηb,c, Pb > 0 1 ηb,d, Pb ≤ 0 (3.24) where ηb,c and ηb,d1 are battery efficiency during charging and discharging mode,
respec-tively.
As was previously mentioned in the section 2.4, the equation 3.24 cannot be exactly handled in the framework of the thesis because it is not continuously differentiable at zero. Moreover, in reality, the battery efficiency increases when the power flow is low [39]. Consequently, the following continuous battery efficiency model with respect to battery flow is used to approximate efficiency of battery in the thesis.
ηb = ηb,c +ηb,d1 2 + ηb,c− ηb,d1 π arctan 100 Pb Pb,max (3.25) where 100 Pb
Pb,max means when the battery power is maximum, the efficiency is 99.3% of the
nominal efficiency. Figure 3.5 shows how equation 3.25 approximates equation 3.24 (The scale is not correct).
Although the battery terminal voltage is considered constant in the thesis, but it is worth mentioning some previous works in the case of further development of the battery module. Battery terminal voltage varies as a function of current, capacity, state-of-charge, and other factors requiring a dynamic model to characterize the voltage at a given time.
Figure 3.5: Model of battery efficiency respect to battery power flow
Different models have been utilized in literature. Guasch et al proposed a dynamic battery model for photovoltaic applications by using automatic parameter extraction techniques and the solving of numerical calculation problems [40]. Manwell et al described kinetic battery model especially for Lead Acid which is used by HOMER energy [41]. Finally, Tremblay et al modeled the battery based on a voltage model in which the discharge curve of a cell is divided into three main zones, including exponential, nominal and over-discharge zone [42]. In fact, all of these models can be implemented in the framework but they are not used since the focus of the thesis is on high-level optimization of microgrid.
3.1.5
AC/DC converters
The inverter in a system converts DC power to AC power and rectifier changes in the opposite direction. The efficiency of these conversions is the most important part of evaluating the system performance. Manufacturers’ datasheets contain a variety of infor-mation essential to the successful application of a PV inverter, including AC voltage, AC frequency, maximum AC power and current, acceptable DC voltage range, maximum DC power and current, DC start up voltage etc.
However, the inverter’s “power conversion performance” or efficiency is often provided as a single peak efficiency value, which can be misleading, and sometimes as a “California
3.1. Electrical system Energy Commission (CEC)” or “European” weighted efficiency value. In low-level system modeling, this single value is not sufficient in order to analyze a power flows and a more accurate model should be utilized. For instance, Sandia inverter model is an implementa-tion of the empirical model described in King 2007 [43] which estimated AC power output of the inverter based on DC power input, DC terminal voltage and parameters of inverter. Since the main focus is one optimization of high-level system, the DC components (PV and battery) terminal voltages are constant and the Sandia inverter model is reduced to a relation between power input and output. For the sake of simplicity, the inverter’s power conversion performance or single efficiency is used to convert DC power generated from PV panels to AC and the battery converter system is assumed that included inside battery system. Hence,
Pinv = ηinvPp (3.26)
where ηinv is inverter efficiency and Pinv is AC power out put of inverter.
3.1.6
Diesel generator
A diesel engine is a component which burns fuel to produce electricity based on a variable efficiency at different loads. The method used for calculation of this efficiency is the same in HOMER [44]. The efficiency of diesel generator, ηgen, is the produced electricity, Pgen,
divided by the chemical energy consumed by the plant. ηgen=
Pgen
˙
mgenLHVgen
(3.27) where ˙mgen is the mass flow rate of the fuel and LHVgen is the lower heating value of
the fuel. ˙mgen is related to rated capacity of the generator, Ygen and the power output.
Homer uses a linear function as below, ˙
mgen= F0Ygen+ F1Pgen (3.28)
where F0 and F1 are diesel consumption coefficients.
3.1.7
Electrical Network
The electrical grid sets the voltage of the AC part of the system and fulfills the customer electrical demand whenever is required. The net power of electrical grid is important since
it determines the associated operation cost of the system. In the case of different tariffs for selling and purchasing electricity, the net power of electrical grid needs to be split into the positive Pnet,pos and negative part Pnet,neg in a differentiable way. Pnet,pos and Pnet,neg
are power consumed from grid and injected to the grid, respectively. Hence,
Pnet = Pnet,pos+ Pnet,neg (3.29)
Since the framework can handle constraints, two attributes, min = 0 and max = 0, are added to Pnet,pos and Pnet,neg, respectively. Then instead of calculating Pnet,neg, it is used
as a input variable and a quadratic penalty (RP2
net,neg) is added to the objective function.
When Pnet,neg is positive, the minimum value of the penalty is zero and based on equation
3.29, the Pnet,pos is equal to Pnet. In the case of negative value for Pnet, the minimum
value of the penalty is at the point where Pnet,pos is zero, then Pnetis Pnet,neg. By choosing
infinitesimal value for R, the effect of this penalty in the objective function becomes negligible.
3.2
Thermal system
In this thesis, the thermal supply is the amount of hot water demanded by customer which can be provided either by thermal network or cogeneration plant. Thermal components used in this work, are obtained from one of Modelon libraries except for the heat exchanger and thermal tank. In this library, for the sake of simplicity, pressure is not modelled; as a result, for new additional components, the pressure modeling is ignored. In the following section, each thermal component model is explained.
3.2.1
Thermal tank
To draw a simple analogy with electrical system, the combination of a tank and its con-nected heat exchanger works as the battery in the electrical system. Hence, thermal tank is a system used to increase the dynamics of the thermal system which can store and distribute heat energy based on the microgrid control decision.
For the sake of simplicity, uniform temperature is assumed for the tank being affected by transferring heat flow with heat exchanger and thermal loss due to temperature difference
3.2. Thermal system with surrounding. Hence,
ρVtCp
dTt
dt = Qt,he− Kt(Tt− Tamb) (3.30) where Vtis the volume of the tank, ρ and Cp are water density and specific heat capacity,
respectively, Ttis tank internal temperature, Qt,heis the heat transfer with heat exchanger
and Kt is the heat loss coefficient.
3.2.2
Heat exchanger
A heat exchanger is a system for transferring heat between two or more fluids. In this thesis, the heat exchanger is used to transfer energy to the thermal tank for later con-sumption. To model the heat exchanger, heat exchanger effectiveness method, number of transfer units (NTU), is used since it is straightforward when more than one of the inlet and outlet temperaturse of the heat exchanger is unknown [45]. It should be noted due to the assumption of uniform temperature in the tank, there is only one unknown temperature in the system; as a result, logarithmic mean temperature difference (LMTD) can be used as well but to open the possibility of further development, the NTU method is chosen.
Based on the NTU method, the number of transfer units (NT U) is indicative of the size of the heat exchanger and is calculated as below,
N T U = UheAhe Cmin,he
(3.31) where Uhe is heat transfer coefficient, Ahe is heat transfer area and Cmin,heis the minimum
of the product of mass flow rate to specific heat capacity between two flows in heat exchanger. The uniform tank temperature assumptions means the heat exchanger tubes are inside a batch of fluid in the tank. Therefore, it is assumed the capacity flow rate on the tank side is far larger than on the heat exchanger side. Based on this assumption, Cmin,he is always for the heat exchanger side.
Cmin,he = ˙mheCp (3.32)
In this method, the heat exchanger effectiveness he is defined as the ratio between the
on uniform tank temperature it can be written as below, he = Qt,he Qmax,he = Tin,he− Tout,he Tin,he− Tt (3.33) where Tin,he and Tout,he are inlet and outlet temperature respectively. By considering the
assumption for Cmin,he, the relation between he and NT U is in the form of[45],
he = 1 − exp −N T U (3.34)
3.2.3
Pumps and volumes
Since pressure is not modeled in this work, the only advantages to use pumps is to put limitation for the mass flow rate between different components in optimization problem (section 2.3.1).
Volumes represent the volume of water in pipes and the distribution between different lines. They do not store water meaning that the sum of all incoming and the sum of all outgoing flows are equal. For instance, in a volume with two inputs (subscript by 1 and 2) and one output (subscript by 3), the mass flow continuity equation is in the form of,
˙
m3 = ˙m1+ ˙m2 (3.35)
Moreover, volumes can be dynamic (3.36a) or static (3.36b) meaning that the internal temperature can be affected or not by incoming and outgoing enthalpy (h) of the flows.
ρVvCp
dTv
dt = ˙m1h1+ ˙m2h2− ˙m3h3 (3.36a) 0 = ˙m1h1+ ˙m2h2− ˙m3h3 (3.36b)
3.2.4
Cogeneration plant
In the thesis, since the focus is on high-level optimization, the details of the power pro-duction are not considered. The cogeneration plant is assumed to produce both electrical power and thermal heat based on variable efficiency with respect to the level of its load. Lcog denotes the load of the cogeneration plant which is 100 in full-load and 0 in no-load
3.3. Modeling in Dymola f2). Hence, ηtot,cog = Qtot,cog− Ql,cog Qtot,cog = f1(Lcog) (3.37a) ηel,cog = Pcog Qt,cog− Ql,cog = f2(Lcog) (3.37b)
Qtot,cog− Ql,cog = Pcog+ Qcog (3.37c)
where Qtot,cogis total heat power from fuel to cogeneration plant at each timestep, Ql,cog is
power loss,Pcog is electrical power of the unit and Qcog is useful thermal power for heating
which transfers energy to water flow as below.
Qcog = ˙mcog(hout,cog− hin,cog) (3.38)
3.3
Modeling in Dymola
In the previous sections, the underling physical theory for all components including, ther-mal and electrical, are explained. As it is elaborated in section 2.1.1, the physical equations of each component are only a part of equation section of a Modelica model. Hence, a complete Dymola model needs connection ports or interfaces. Table 3.1 shows different type of ports used in this thesis. Interfaces define physical properties in the boundaries of each model and connected components have same boundaries by using connectors. Table 3.2 and 3.3 represent the summary of electrical and thermal components, respectively.
Table 3.1: Port definition
Port Name Boundaries
DC voltage, current
AC voltage, current
Heat heat power, temperature Flow mass flow rate, enthalpy
Table 3.2: Electrical components summery
Name Icon Theory Input Output Connectionport type
Weather Section3.1.1 Weatherdata
ψs, θs, DHI, DN I, Tamb -PV tracking Section3.1.2 ψs, θs, DHI, DN I Ig -PV Section3.1.3 Ig,Tamb - DC Battery
pack Section3.1.4 - - AC*
Inverter Section3.1.5 - - DC, AC
Diesel
generator Section3.1.6 - - AC
Electrical
grid Section3.1.7 - - AC