An adaptive discretization MINLP algorithm for optimal synthesis of decentralized energy supply systems

Accepted Manuscript Title: An Adaptive Discretization MINLP Algorithm for Optimal Synthesis of Decentralized Energy Supply Systems Author: Sebastian Goderbauer Bj¨orn Bahl Philip Voll Marco E. Lubbecke ¨ Andr´e Bardow Arie M.C.A. Koster PII: DOI: Reference:

S0098-1354(16)30286-1 http://dx.doi.org/doi:10.1016/j.compchemeng.2016.09.008 CACE 5557

To appear in:

Computers and Chemical Engineering

Received date: Revised date: Accepted date:

14-6-2016 6-9-2016 8-9-2016

Please cite this article as: Sebastian Goderbauer, Bjddotorn Bahl, Philip Voll, Marco E. Lddotubbecke, Andr´e Bardow, Arie M.C.A. Koster, An Adaptive Discretization MINLP Algorithm for Optimal Synthesis of Decentralized Energy Supply Systems, (2016), http://dx.doi.org/10.1016/j.compchemeng.2016.09.008 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

An Adaptive Discretization MINLP Algorithm for Optimal Synthesis of Decentralized Energy Supply Systems Sebastian Goderbauer, Björn Bahl, Philip Voll, Marco Lübbecke, André Bardow, Arie M.C.A. Koster Preprint submitted to Computers & Chemical Engineering

ip t

Highlights

cr

+ A solution algorithm for a nonconvex nonlinear synthesis problem is proposed.

us

+ We use an iterative interaction between approximate MIPs and decomposable NLPs.

an

+ An iteratively adapted discretization forms the basis of the approximate problem.

M

+ The algorithm outperforms MINLP solvers and pure linearization approaches.

Ac ce pt e

d

+ Real-world-based test instances used in the computational study are online available.

Page 1 of 39

ip t

An Adaptive Discretization MINLP Algorithm for Optimal Synthesis of Decentralized Energy Supply SystemsI

cr

Sebastian Goderbauera,b , Bj¨orn Bahlc , Philip Vollc , Marco E. L¨ ubbeckeb , Andr´e Bardowc , Arie M.C.A. Kostera,∗ a Lehrstuhl

II f¨ ur Mathematik, RWTH Aachen University, 52056 Aachen, Germany Research, RWTH Aachen University, 52072 Aachen, Germany c Chair of Technical Thermodynamics, RWTH Aachen University, 52056 Aachen, Germany

an

us

b Operations

Abstract

Decentralized energy supply systems (DESS) are highly integrated and complex

M

systems designed to meet time-varying energy demands, e.g., heating, cooling, and electricity. The synthesis problem of DESS addresses combining various types of energy conversion units, choosing their sizing and operations to max-

d

imize an objective function, e.g., the net present value. In practice, invest-

te

ment costs and part-load performances are nonlinear. Thus, this optimization problem can be modeled as a nonconvex mixed-integer nonlinear program-

Ac ce p

ming (MINLP) problem. We present an adaptive discretization algorithm to solve such synthesis problems containing an iterative interaction between mixedinteger linear programs (MIPs) and nonlinear programs (NLPs). The proposed algorithm outperformes state-of-the-art MINLP solvers as well as linearization approaches with regard to solution quality and computation times on a test set obtained from real industrial data, which we made available online. Keywords: Mixed-Integer Nonlinear Programming; Decentralized Energy

Supply System; Synthesis; Structural Optimization; Adaptive Discretization

I Funded

by the Excellence Initiative of the German federal and state governments. author Email address: [email protected] (Arie M.C.A. Koster)

∗ Corresponding

Preprint submitted to Computers & Chemical Engineering

September 6, 2016

Page 2 of 39

1. Introduction

2

We propose an adaptive discretization algorithm for the superstructure-based

3

synthesis of decentralized energy supply systems (DESS). The proposed opti-

4

mization-based algorithm employs discretization of the continuous decision vari-

5

ables. The discretization is iteratively adapted and used to obtain valid noncon-

6

vex mixed-integer nonlinear program (MINLP) solutions within short solution

7

time.

us

cr

ip t

1

8

DESS can consist of several energy conversion components (e.g., boilers and

10

chillers) providing different utilities (e.g., heating, cooling, electricity). DESS

11

are highly integrated and complex systems due to the integration of different

12

forms of energy and their connection to the gas and electricity market as well as

13

to the energy consumers. The application of DESS encompasses, e.g., chemical

14

parks (Mar´echal & Kalitventzeff, 2003), urban districts (Mar´echal et al., 2008;

15

Jennings et al., 2014) and building complexes (Arcuri et al., 2007; Lozano et al.,

16

2009). Energy costs usually match the companies’ profits in magnitude (Drumm

17

et al., 2013). Thus, optimally designed decentralized energy supply systems can

18

lead to a considerable increase of profits.

M

d

te

Ac ce p

19

an

9

20

The target of optimal synthesis of DESS is the identification of an (economi-

21

cally) optimal structure (which types of equipment and how many units?) and

22

optimal sizing (how big?), while simultaneously considering the optimal op-

23

eration of the selected components (which components are operated at which

24

level at what time?) (Frangopoulos et al., 2002). These three decision levels

25

could be considered sequentially. However, the levels influence each other, thus

26

only a simultaneous optimization will find a global optimal solution. In this

27

paper, we consider the simultaneous optimization using superstructure-based

28

synthesis. A superstructure needs to be predefined and consists of a superset of

29

possible components, which can be selected within the synthesis of the DESS. If

30

the superstructure is chosen too small, optimal solutions could be excluded, if

2

Page 3 of 39

the superstructure is chosen too large, computational effort become prohibitive.

32

Therefore some of the authors proposed a successive superstructure expansion

33

algorithm (Voll et al., 2013b).

ip t

31

34

The synthesis of DESS contains binary decisions for the selection of energy

36

conversion components as well as the on/off status in the operation of each

37

component. Combined with nonlinear part-load performance of the energy con-

38

version components, nonlinear economy-of-scale effects in the investment cost

39

curves and strict energy balances, the synthesis of DESS leads in general to

40

a nonconvex MINLP (Bruno et al., 1998). Typically, an economic objective

41

function is considered, e.g., the net present value is maximized or the total an-

42

nualized costs are minimized, furthermore also ecologic objective functions can

43

be considered (Østergaard, 2009).

M

44

an

us

cr

35

Metaheuristic optimization approaches have been proposed for the synthesis of

46

DESS: Evolutionary algorithms were proposed for superstructure-free linearized

47

synthesis as well as superstructure-based MINLP synthesis (Dimopoulos & Fran-

48

gopoulos, 2008; Voll et al., 2012). Stojiljkovic et al. (2014) proposed a heuristic

49

for structural decisions and solved an mixed-integer linear program (MILP) for

Ac ce p

te

d

45

50

operation decision. These heuristic approaches do not provide any measure of

51

optimality.

52

53

To allow rigorous optimization, mostly linearized approaches are considered for

54

synthesis of practically relevant problems. In the resulting MILPs, the non-

55

linearities are approximated by piecewise-linearized functions. First, Papoulias

56

& Grossmann (1983) linearized the investment cost functions, the nonlinear

57

operation conditions are modeled as discrete, but fixed operation conditions.

58

Continuous operation decision with constant efficiency is addressed by Lozano

59

et al. (2009) for MILP synthesis of energy supply systems in the building sec-

60

tor using fixed capacities. Voll et al. (2013b) proposed an MILP model ac-

61

counting for piecewise-linearized part-load dependent operation conditions and 3

Page 4 of 39

piecewise-linearized investment costs for continuous component sizing. Recently,

63

Yokoyama et al. (2015) modeled the structure decision with integer variables for

64

the type and discrete sizes of components, thus, modeling the nonlinear invest-

65

ment cost curve is not required. The operation power is modeled as linear

66

function within allowed operation ranges.

cr

ip t

62

67

The solution of the linarization approaches only results in approximated solu-

69

tions. However solving the MINLP of superstructure-based synthesis is compu-

70

tationally demanding. First, an MINLP model for the operation of DESS was

71

considered by Prokopakis & Maroulis (1996). The model takes into account

72

the nonlinear size- and load-dependent components performance. Papalexandri

73

et al. (1998) and Bruno et al. (1998) generalized the MINLP formulation to

74

the optimal synthesis of DESS. Due to the complexity of the problem, only one

75

component of each type is considered in the superstructure and the demand

76

is considered by a single load case. An MINLP model considering multiple,

77

detailed components as well as multiple load cases for the demand profile has

78

been proposed by Varbanov et al. (2004, 2005). To solve the resulting large

79

MINLP, nonlinearities of part-load performance are predefined in an iterative

80

loop and internally MILPs are solved. Chen & Lin (2011) solved an MINLP

Ac ce p

te

d

M

an

us

68

81

for a steam-generation plant, the nonlinearities of part-load performance are

82

optimized, nevertheless the model considers steam as a single demand type.

83

The problem of integrated optimization of DESS and process system commonly

84

results in large-scale MINLPs. Recently Zhao et al. (2015) decomposed the

85

integrated MINLP of optimal operation of DESS and process system into an

86

MILP and NLP problem and the variables are exchanged between both prob-

87

lems. Moreover, Tong et al. (2015) proposed a discretization approach for the

88

MINLP of optimal operation of DESS and process system. Further discretiza-

89

tion approaches for solving nonconvex MINLP problems with different practical

90

applications are discussed in Section 3

91

92

In this paper, MINLP solutions are obtained by an adaptive discretization al4

Page 5 of 39

gorithm for the nonlinear synthesis problem of DESS. (Commercial) MINLP

94

solvers such as BARON (Tawarmalani & Sahinidis, 2005) reach computational

95

limits for relative small test cases of the considered MINLP, accounting for

96

nonlinear investment cost and multivariate nonlinear part-load dependent op-

97

eration performance. We developed a problem-tailored adaptive discretization

98

algorithm to obtain valid solutions of the MINLP within short solution time.

99

The algorithm discretizes the continuous component size within bounds given

100

by practically available component size limits. The whole range of size can

101

be selected for each type of component, since the discretization is iteratively

102

adapted. Thus, the algorithm does not require predefining discrete sizes of the

103

components in the superstructure. Moreover, the operation of each component

104

for each load case is discretized with finer steps depending on the part-load per-

105

formence of each type of energy conversion component. Thus, various energy

106

conversion components with different capacities and with corresponding invest-

107

ment and maintenance costs can be selected and adjusted to meet the energy

108

demands in each load case.

109

We state our MINLP model of the DESS in Section 2. In Section 3, we describe

110

the proposed adaptive discretization algorithm. In Section 4, we apply the al-

111

gorithm to a test set of a real-world example. Solutions and performance are

Ac ce p

te

d

M

an

us

cr

ip t

93

112

compared to a standard MINLP solver as well as state-of-the-art linearization

113

approaches with MILP models.

114

2. Optimization Models for Decentralized Energy Supply Systems

115

In this section, we present an MINLP model for optimal synthesis of DESS

116

(Section 2.2) as well as a piecewise-linearized model (Section 2.3), which we use

117

as benchmark for our adaptive discretization algorithm. First of all, in Section

118

2.1, notations of parameters, decisions, and the optimization problem as a whole

119

are given.

5

Page 6 of 39

Boiler

CHP engine

Legend:

Electricity Heat

Power supply

Electrictiy demands

. Cooling E cool demands

. el

E

us

. heat Heat E demands

Turbo chiller

cr

Absorption chiller

ip t

Cooling

Natural gas hook-up

Figure 1: Example of a decentralized energy supply system with exactly one unit of each

an

considered type of equipment.

2.1. Equipment, Parameters, and Decisions

121

The set of energy conversion units, which can be set up to meet the demands, is

122

denoted by superstructure S = B ∪˙ C ∪˙ T ∪˙ A and encompasses a set of boilers

123

B, a set of combined heat and power engines C, a set of turbo-driven compressor

124

chillers T and a set of absorption chillers A (Figure 1). Further equipment

125

could be included, but we focus here on the problem introduced in our earlier

126

work Voll et al. (2013b). All units s ∈ S in the superstructure are not further

127

specified than their type of equipment. Note, that an optimal DESS is likely to

Ac ce p

te

d

M

120

128

contain multiple units of one type which is in strong contrast to classical process

129

synthesis problems (Farkas et al., 2005).

130

The set of load cases considered for the operation of the DESS is denoted by L.

131

The length of load case ℓ ∈ L is denoted by ∆ℓ ≥ 0. Furthermore, E˙ ℓheat ≥ 0,

132

E˙ ℓcool ≥ 0, and E˙ ℓel ≥ 0 denote the demands of heating, cooling, and electricity,

133

which have to be satisfied with equality by the DESS in every load case ℓ ∈ L.

134

For each unit s ∈ S, its continuous size V˙ sN has to be determined. The size

135

V˙ sN specifies the maximum (nominal) output energy and has to be between a

136

minimum size V˙ sN,min and a maximum size V˙ sN,max . For combined heat and

137

power (CHP) engines, the output is not unique (heat and electricity). In this

138

case, the size refers to the maximum heat output. The investment cost of unit

6

Page 7 of 39

s ∈ S depends on its size V˙ sN and is given by the nonlinear function Is (V˙ sN ).

140

Further, maintenance costs are considered as constant factors ms in terms of

141

investment costs.

142

The output power of unit s ∈ S at load case ℓ ∈ L is to be determined and is

143

denoted by V˙ sℓ . Again, for CHP, the output power refers to the heat output.

144

The nonlinear function V˙ sℓel (V˙ sℓ , V˙ sN ) describes the electricity output of a CHP

145

s ∈ C ⊆ S. For each unit s ∈ S operated in load case ℓ ∈ L, a minimum

146

part-load operation is required. Thus, the condition αsmin V˙ sN ≤ V˙ sℓ ≤ V˙ sN with

147

minimum part-load factor 0 ≤ αsmin ≤ 1 has to hold. If s ∈ S is not operated in

148

load case ℓ ∈ L, we set V˙ sℓ = 0. The input needed to generate the output V˙ sℓ is

149

described by the nonlinear part-load performance function U˙ s (V˙ sℓ , V˙ sN ).

an

us

cr

ip t

139

Parameters pgas,buy , pel,buy , and pel,sell denote the purchase price of gas and

M

electricity, and the selling price of electricity from and to the grid. To compute the objective value of a feasible DESS, i.e., the net present value, the parameter

d

APVF(i, γ CF ) :=

CF

(i + 1)γ − 1 i · (i + 1)γ CF

denotes the present value factor and depends on discount rate i and cash flow

151

time γ CF .

Ac ce p

152

te

150

153

The equipment models including the analytical equations of the unit’s input-

154

output relations and all parameters can be found in Appendix A.

155

2.2. MINLP Formulation

156

Variables. For every unit s ∈ S, the variable ys ∈ {0, 1} denotes whether the

157

unit is chosen and the continuous variable V˙ sN ≥ 0 denotes the (nominal) size of

158

unit s. The variable δsℓ ∈ {0, 1} denotes the on/off-status and the continuous

159

variable V˙ sℓ ≥ 0 denotes the output power of unit s ∈ S in load case ℓ ∈ L.

160

Furthermore, the continuous variables U˙ ℓel,buy ≥ 0 and V˙ ℓel,sell ≥ 0 denote the

161

bought and sold electricity power from and to the grid in load case ℓ ∈ L.

162

Formulation. A mixed-integer nonlinear programming formulation for the con-

163

sidered problem for optimal synthesis of DESS is given by (1) – (12). 7

Page 8 of 39

( ∆ℓ · pel,sell · V˙ ℓel,sell − pel,buy · U˙ ℓel,buy

ℓ∈L

∑

− pgas,buy · −

∑

s∈B∪C

V˙ sℓ = E˙ ℓheat +

s∈B∪C

∑

s∈S

∑

δsℓ · U˙ s (V˙ sℓ , V˙ sN )

s∈A

∀ℓ∈L

an

∑

)

(1)

] ∑ N ˙ ms · Is (Vs ) · ys − Is (V˙ sN ) · ys

s∈S

s.t.

δsℓ · U˙ s (V˙ sℓ , V˙ sN )

cr

[∑

us

max APVF(i, γ CF ) ·

ip t

164

∀ℓ∈L

V˙ sℓ = E˙ ℓcool

s∈A∪T

(2) (3)

M

∑ ∑ U˙ ℓel,buy + δsℓ · V˙ sel (V˙ sℓ , V˙ sN ) = E˙ ℓel + δsℓ · U˙ s (V˙ sℓ , V˙ sN ) + V˙ ℓel,sell ∀ℓ∈L

(4)

∀s∈S

(5)

∀ s ∈ S, ℓ ∈ L

(6)

∀ s ∈ S, ℓ ∈ L

(7)

V˙ sℓ ≥ αsmin · V˙ sN − (1 − δsℓ ) · αsmin · V˙ sN,max

∀ s ∈ S, ℓ ∈ L

(8)

ys ≥ δsℓ

∀ s ∈ S, ℓ ∈ L

(9)

s∈C

s∈T

Ac ce p

V˙ sℓ ≤ V˙ sN

te

V˙ sℓ ≤ δsℓ · V˙ sN,max

d

V˙ sN,min ≤ V˙ sN ≤ V˙ sN,max

ys ∈ {0, 1}, V˙ sN ≥ 0

∀s∈S

(10)

δsℓ ∈ {0, 1}, V˙ sℓ ≥ 0

∀ s ∈ S, ℓ ∈ L

(11)

U˙ ℓel,buy , V˙ ℓel,sell ≥ 0

∀ℓ∈L

(12)

165

Objective. In objective (1), the net present value NPV := APVF(i, γ CF ) ·

166

RCF − I is maximized. The NPV is calculated from the present value factor

167

APVF(i, γ CF ), the net cash flow RCF and the total investments I. The net cash

168

flow RCF are the annual revenues from sold electricity V˙ ℓel,sell minus the cost

169

for electricity U˙ ℓel,buy bought from the grid and secondary energy U˙ s (V˙ sℓ , V˙ sN )

170

consumed by the boilers and CHP engines as well as maintenance costs. 8

Page 9 of 39

Constraints. Constraints (2) – (4) ensure that the demands for heating Eℓheat ,

172

cooling Eℓcool and electricity Eℓel are fulfilled with equality in every load case

173

ℓ ∈ L by the DESS. Constraints (5) restrict the size VsN to be in the tech-

174

nically allowed range [VsN,min ; VsN,max ]. Constraints (6) – (8) force Vsℓ = 0,

175

if δsℓ = 0 and, otherwise, limit Vsℓ to the operation interval [αsmin · VsN ; VsN ].

176

Constraints (9) ensure that a unit is chosen, if it is used in at least one load case.

cr

ip t

171

us

177

We note that the formulation is nonlinear due to the equipment models of the

179

units (Appendix A) and bilinear terms in (1), (2), (4) as well as nonconvex due

180

to the investment cost function Is (V˙ sN ) (Appendix A) and nonlinear equality

181

constraints (2) and (4).

182

2.3. Benchmarking to Piecewise-linearized Approach

183

The MINLP synthesis model (1) – (12) is commonly linearized for practical ap-

184

plications (Section 1). The solution obtained by the approximated MILP is in

185

general not feasible for the nonlinear model (Bruno et al., 1998) (Section 4.2). In

186

this section, we present an approach to obtain feasible solutions of the MINLP

187

based on a solutions of the MILP. The feasible MINLP solution based on the

188

MILP result is considered as benchmark for our adaptive discretization algo-

te

d

M

an

178

rithm (Section 3). Since, as explained above, a one-to-one comparison between

190

MINLP and MILP solutions is not possible, we think that the presented analy-

191

sis provides an insightful comparison between previous work and the algorithm

192

proposed in this work.

193

The MILP stated by Voll et al. (2013b) with piecewise-linearized functions for

194

the nonlinear investment cost curves and piecewise-linearized part-load opera-

195

tion curves are employed to compute the MILP solution. To obtain a feasible

196

∗ ˙∗ MINLP solution based on the solution δsℓ , Vsℓ , ys∗ , V˙ sN∗ of the MILP, we solve

197

MINLPlin,feas (13) – (15).

Ac ce p 189

9

Page 10 of 39

    V˙ sN∗ − V˙ sN V˙ sℓ∗ − V˙ sℓ ∑  +   min |L| · V˙ sN∗ V˙ sℓ∗ s∈S s∈S,ℓ∈L: ∗ δsℓ =1

s.t. (2) − (12)

(13)

ip t

∑

∗ ys = ys∗ , δsℓ = δsℓ

cr

(14)

∀s ∈ S, ℓ ∈ L

(15)

∗ The selected structure of the DESS ys∗ and the on/off status of the equipment δsℓ

199

defined by the MILP is kept fixed. The objective function reflects minimizing

200

the difference between the solution values of the MILP and the MINLP, in

201

the solution space of the MINLP. Thus, feasible solutions for the MINLP are

202

obtained which are ‘near’ the given solution of the MILP. The difference measure

203

is defined by the sum of the normalized differences in optimal design V˙ sN and

204

operation V˙ sℓ .

205

3. Adaptive Discretization Algorithm

206

Solving the nonconvex MINLP (1) – (12) with state-of-the-art solvers like

207

BARON (Tawarmalani & Sahinidis, 2005) leads to unsatisfying results. For

208

several nontrivial test instances, it is even hard for solvers to compute a feasi-

Ac ce p

te

d

M

an

us

198

209

ble solution (Section 4). Thus, the need of a problem-specific solution method

210

providing primal MINLP solutions is evident.

211

212

It is a common approach to discretize (continuous) variables in a nonconvex

213

nonlinear program to approximate it with an easier to solve mixed-integer lin-

214

ear one (e.g., Leyffer et al. (2008), Pham et al. (2009), Geißler et al. (2011),

215

Gupte et al. (2013), Kolodziej et al. (2013), Yue & You (2014)). As an obtained

216

solution might not be feasible for the original MINLP, we extend the discretiza-

217

tion approach. Given an approximate solution, we fix selected solution-specific

218

decisions (i.e., unit sizes and on/off statuses in the load cases), and solve the

219

remaining NLP using a decomposition to arrive at a primal MINLP solution. To

220

have a computational tractable approximation, only a few discretization points 10

Page 11 of 39

for the unit size and operation are used. However, to ensure a certain accuracy

222

in our discretization algorithm, this two step algorithm is embedded in a loop

223

of refinements of the discretization. At every iteration, the discretization grid

224

is contracted and shifted in the direction of the discrete unit size chosen in the

225

previous iteration, keeping the size of the discretized MILP.

cr

ip t

221

226

The discretized problem formulation of the nonconvex MINLP (1) – (12) is

228

described in Section 3.1, followed by the procedure to form a feasible MINLP

229

solution using the approximate solution in Section 3.2. The adaptive part of the

230

discretization algorithm is specified in Section 3.3. In the end and putting every-

231

thing together, Section 3.4 provides a description of the adaptive discretization

232

algorithm as a whole and some further comments.

233

3.1. Discretized Problem

234

To develop an approximation via discretization, we discretize all continuous

235

variables with nonlinear dependencies. In MINLP (1) – (12) this involves the

236

unit’s size V˙ sN ≥ 0 and operation V˙ sℓ ≥ 0.

te

d

M

an

us

227

Unit size. If unit s ∈ S is chosen, i.e., ys = 1 holds, we have to choose a size V˙ sN

Ac ce p

in the interval [V˙ sN,min , V˙ sN,max ]. The investment cost function Is (V˙ sN ) depends on the unit’s size and is nonlinear (Appendix A). We discretize the range of the continuous variable V˙ sN ∈ [V˙ sN,min , V˙ sN,max ] by dividing the interval into ksmax (equidistant) discrete sizes

N,val N,val N,val V˙ s1 < V˙ s2 < . . . < V˙ sk max s

237

N,val ˙ N,val , Vskmax ≤ V˙ sN,max and ksmax ∈ N an odd number. These with V˙ sN,min ≤ V˙ s1

238

N,val N V˙ sk ∈ {0, 1} denote whether are parameters and the related variables V˙ sk

239

N,val unit s ∈ S is set up with the k-th discrete size V˙ sk . Thus, we transform every

240

N continuous variable V˙ sN in several binary variables V˙ sk . This implies that we

241

do not need the nonlinear investment cost function in the discretized problem

242

N,val anymore, because for each discrete size V˙ sk we can compute its investment

243

N,val val costs Isk := Is (V˙ sk ) in advance.

s

11

Page 12 of 39

Unit operation. Together with the unit’s size, one can compute the input energy needed to get a desired operation output V˙ sℓ using the nonlinear part-load

ip t

N,val performance functions U˙ s (V˙ sℓ , V˙ sN ) (Appendix A). If V˙ sk is chosen out of the

discrete unit sizes and the unit is switched on, the possible output V˙ sℓ of this

cr

N,val N,val unit is bound by the size V˙ sk and a minimal possible part-load αsmin V˙ sk

with 0 ≤ αsmin ≤ 1. Again, we discretize the range of the continuous vari-

us

N,val ˙ N,val max able V˙ sℓ ∈ [αsmin V˙ sk , Vsk ] by dividing the interval into jskℓ (equidistant)

discrete operations

an

N,val val val val ˙ N,val αs · V˙ sk =: V˙ skℓ1 < V˙ skℓ2 < . . . < V˙ skℓj max := V sk skℓ

max max with jskℓ ∈ N and jskℓ ≥ 2. The variable V˙ skℓj ∈ {0, 1} denotes whether unit N,val val s ∈ S with size V˙ sk has the j-th discrete operational output V˙ skℓj in load

M

case ℓ ∈ L. We name

N,val N,val el,val val val val , V˙ sk , V˙ sk := U˙ s (V˙ skℓj ) and V˙ skℓj ) U˙ skℓj := Vsel (V˙ skℓj

the values of the employed part-load performance function at the corresponding

245

discrete size and discrete operation.

246

The discretization grid of size and operation and its notation is summarized in

247

Figure 2.

Ac ce p

te

d

244

Since fulfilling the energy demands in heating, cooling, and electricity with equality is a requirement of our problem (constraints (2) – (4)), we want to enable equality in the energy balances of our discretized approximation as well (cf. Remark 1 in Section 3.4). A pure discretization with binary variables

V˙ skℓj does not allow this in general. Therefore, we piecewise linearize the part-

val load performance functions U˙ s and V˙ sel using V˙ skℓj as supporting points. We

introducing a continuous variable cont V˙ skℓj ≥0

with

val,diff cont val val V˙ skℓj ≤ V˙ skℓj+1 − V˙ skℓj =: V˙ skℓj

and parameters lin U˙ skℓj :=

val val U˙ skℓj+1 − U˙ skℓj V˙ val,diff

,

skℓj

el,lin V˙ skℓj :=

el,val el,val V˙ skℓj+1 − V˙ skℓj V˙ val,diff skℓj

12

Page 13 of 39

val val for each simplex, i.e., line segment between V˙ skℓj and V˙ skℓj+1 . Note that in the

249

proposed approach the discretization of the unit size remains a pure one, where

250

we do not add further continuous variables or piecewise linearize anything there.

ip t

248

251

Using the specified discretization and linearization, we are able to approximate

253

the original MINLP (1) – (12) with the following mixed-integer linear program

254

(16) – (26) (discretized MIP). For better readability, the limits of the indices

255

max k = 1, 2, . . . , ksmax and j = 1, 2, . . . , jskℓ are not mentioned explicitly in the

256

following formulation and sections. Indices s and ℓ without any set information

257

mean s ∈ S and ℓ ∈ L.

an

us

cr

252

k=2

k=4

···

···

···

Ac ce p

···

te

N V˙ sk

k=3

VsN,max

d

k=1

M

VsN,min .

N αmin V˙ sk s

ℓ=1 ℓ=2

··· j=2

j=1 ···

Figure 2: Discretization and notation of unit size decisions (top) and operation decisions (bottom).

13

Page 14 of 39

( ∆ℓ · pel,sell · V˙ ℓel,sell − pel,buy · U˙ ℓel,buy

ℓ∈L

− pgas,buy ·

∑ ∑( )) val lin cont U˙ skℓj · V˙ skℓj + U˙ skℓj · V˙ skℓj

s∈B∪C k,j

−

∑

ms ·

val Isk

] ∑ N val ˙ N ˙ · Vsk − · Vsk Isk

s,k

val cont V˙ skℓj · V˙ skℓj + V˙ skℓj

s∈B∪C k,j = E˙ ℓheat

+

∑∑(

)

val lin cont U˙ skℓj · V˙ skℓj + U˙ skℓl · V˙ skℓj

s∈A k,j

∑ ∑(

) val cont V˙ skℓj · V˙ skℓj + V˙ skℓj = E˙ ℓcool

s∈A∪T k,l

(16)

us

∑ ∑(

s,k

)

an

s.t.

ip t

[∑

cr

max APVF(i, γ CF ) ·

∀ ℓ ∈ L (17) ∀ ℓ ∈ L (18)

M

∑ ∑ ( el,val el,lin ˙ cont ) U˙ ℓel,buy + V˙ skℓj · V˙ skℓj + V˙ skℓj · Vskℓj

s∈C k,j ∑∑( ) val lin cont U˙ skℓj · V˙ skℓj + U˙ skℓj · V˙ skℓj ∀ ℓ ∈ L (19) = E˙ ℓel + V˙ tel,sell +

d

s∈T k,j

te

val,diff ˙ cont V˙ skℓj ≤ V˙ skℓj · Vskℓj ∑ N ≤1 V˙ sk k

∀ s ∈ S (21)

Ac ce p

∑

∀ s, k, ℓ, j (20)

∀ s, k, ℓ (22)

N V˙ skℓj ≤ V˙ sk

j

N V˙ sk ∈ {0, 1}

∀ s, k (23)

V˙ skℓj ∈ {0, 1}

∀ s, k, ℓ, j (24)

cont V˙ skℓj ≥0

∀ s, k, ℓ, j (25)

U˙ ℓel,buy , V˙ ℓel,sell ≥ 0

∀ ℓ ∈ L (26)

258

The binary variables ys and δsℓ in MINLP formulation (1) – (12) are not needed

259

N and V˙ skℓj together with constraints anymore, since the new binary variables V˙ sk

260

(21) and (22) include their role. Of course, the discretized problem (16) – (26)

261

contains a lot more binary decisions, but the nonlinearities and actually the 14

Page 15 of 39

nonconvex nonlinearities are eliminated in that model. The discretized problem

263

(16) – (26) is a mixed-integer linear program. Moreover, some binary variables

264

N V˙ sk and V˙ skℓj can be eliminated by preprocessing (Section 4.2) depending on

265

the demands.

266

3.2. Nonlinear Feasibility

cr

ip t

262

After computing a solution of (16) – (26), post processing is needed to compute

us

a primal feasible solution of the original MINLP (1) – (12). For this post processing, we fix the unit sizes and load case specific on/off statuses given by N∗ ∗ the approximate solution. Therefore, let V˙ sk ∈ {0, 1} and V˙ skℓj ∈ {0, 1} be the

an

values of the related variables of a given discretized problem solution. For each unit s ∈ S and load case ℓ ∈ L parameters ∑

N∗ V˙ sk

(27)

N,val ˙ N∗ V˙ sk · Vsk

(28)

∗ V˙ skℓj

(29)

M

yspar :=

k

V˙ sN,par :=

∑ k

∑

d

par δsℓ :=

te

k,j

are defined. Fixing the variables ys ∈ {0, 1}, V˙ sN ≥ 0 and δsℓ ∈ {0, 1} with

268

these values in MINLP (1) – (12) implies that all binary variables and all con-

Ac ce p

267

269

straints (5) – (9) linking load case are not needed anymore or become variables

270

bounds. As a consequence, the problem is decomposable in independent non-

271

linear programs, one for every load case ℓ ∈ L, named NLPℓ . Since the equality

272

constraints (2) and (4) are still present, every NLPℓ remains nonconvex. How-

273

ever, as the computational results show in section 4, the independent problems

274

are solved quite fast for the considered test instances. It is not guaranteed that

275

NLPℓ provides a feasible solution, since parts of an approximate solution are

276

fixed. Whether NLPℓ provides a feasible solution depends on the form of the

277

piecewise linearized functions and on the fineness of the discretization. It turns

278

out that more discretization points, on the one hand, enlarge the probability to

279

determine a feasible solution in NLPl but, on the other hand, enlarge the com-

280

putation times of solving NLPl . However, the computational results (Section 15

Page 16 of 39

4.2) show that for all test instances considered in this paper every single NLPℓ

282

was feasible.

283

3.3. Adapting Discretization

284

Solving the discretized problem (16) – (26) followed by suited NLPs for nonlin-

285

ear feasibility, a primal solution of the MINLP (1) – (12) is probably computed.

286

This interaction of MIP and NLPs is incorporated into an iteration loop, as it

287

is described as a whole in Section 3.4. At the end of every iteration step, the

288

discretization grid is adapted based on the just computed solution of the dis-

289

cretized problem in that step. Thereby, nearly the entire spectrum of possible

290

unit sizes is enabled and a greater accuracy in the whole algorithm and therefore

291

better primal solutions are achieved. In the rest of this section, the procedure

292

of adapting the discretization is described in detail. Proceed to Section 3.4 for

293

an overview of the whole adaptive discretization algorithm.

M

an

us

cr

ip t

281

294

N,val N,val N,val Let V˙ s1 < V˙ s2 < . . . < V˙ sk max be the discretization of the size of unit s ∈ S

d

s

in a certain iteration step and let 1 ≤ ks∗ ≤ ksmax be the index of the chosen

te

N,val discrete size of unit s ∈ S, i.e., V˙ sN,par = V˙ sk holds (cf. equation (28)). We ∗ s

Ac ce p

are faced with three different cases depending on ks∗ and V˙ sN,par : Case 1 (interior point):

ks∗ ̸∈ {1, ksmax }

Case 2 (interval limit):

N,val ̸∈ {V˙ sN,min , V˙ sN,max } ks∗ ∈ {1, ksmax } and V˙ sk ∗

Case 3 (bound):

ks∗

s

∈

{1, ksmax }

N,val and V˙ sk ∈ {V˙ sN,min , V˙ sN,max } ∗ s

295

Case 1 (interior point). For the case ks∗ ̸∈ {1, ksmax } the chosen discrete point

296

lies not at the end of the discretization interval. We then contract the dis-

297

N,val crete grid points in the direction of the chosen size V˙ sk . More precisely, the ∗

298

N,val N,val equidistant division of the interval [Vsk ∗ −1 , Vsk ∗ +1 ] gets the new and adapted

299

discretization. Notice, since ksmax is odd, V˙ sN,par stays to be a grid point after

300

the adaption.

s

s

s

16

Page 17 of 39

303

N,val N,val V˙ s1 ← Vsk ∗ s −1 ) ( 1 N,val N,val N,val V˙ s2 ← Vsk + Vsk ∗ ∗ −1 s s 2 N,val N,val V˙ s3 ← Vsk∗ s ( ) 1 N,val N,val N,val V˙ s4 ← Vsk ∗ +1 + Vsk ∗ s s 2 V˙ N,val ← V N,val . ∗ s5

ip t

ks∗

cr

302

Example (interior point). For ksmax = 5, ks∗ ̸∈ {1, 5} the adaption effects:

Figure 3: Adaptive Discretization (interior point)

us

301

sks +1

an

Case 2 (interval limit). If ks∗ ∈ {1, ksmax } is the index of one of the end points . N,val of the discretization and V˙ sk ̸∈ {V˙ sN,min , V˙ sN,max } holds, we shift the grid in ∗ s

the direction of the chosen size (and do not contract it). Thus, for ks∗ = 1 (case

s

M

ks∗ = knmax analog) the equidistant division of the interval [ ] N,val N,val N,val ˙ ˙ ˙ ⌉ ⌈ ⌉ ⌈ 2 · Vsk∗ − V ksmax , V ksmax s

2

s

2

results in the adapted discretization (Figure 4a). It is guaranteed that this

305

discretization respects the initial size bounds V˙ sN,min and V˙ sN,max .

306

N,val Case 3 (bound). In the remaining case of ks∗ ∈ {1, ksmax } and V˙ sk ∈ {V˙ sN,min , V˙ sN,max } ∗

te

d

304

s

Here, we assume ks∗ = 1 and

the chosen discrete point lies at a bound.

308

N,val V˙ sk = V˙ sN,min (other case analog). We contract the discretization, respect∗

Ac ce p 307

s

309

ing the bounds V˙ sN,min and V˙ sN,max . The equidistant division of the interval

310

⌈ kmax ⌉ ] is the adapted discretization (Figure 4b). [V˙ sN,min , V˙ N,val s s

2

ks∗ = 1

ks∗ = 1 V˙ sN,min

(a) case 2: chosen (interval limit)

(b) case 3 (bound)

Figure 4: Adaptive Discretization: (a) case 2 and (b) case 3

. 17

.

Page 18 of 39

Example (interval limit, bound). Figure 4 shows examples for the adaption of

312

the discretization in case 2 and 3.

313

3.4. Adaptive Discretization Algorithm

314

The adaptive discretization algorithm is shown in Figure 5 as a whole. The

315

algorithm consists of an iteration loop with a certain stop criterion, for example

316

concerning an iteration or convergence limit. If the criterion is fulfilled, the

317

algorithm terminates and the best MINLP solution found is the output of the

318

algorithm, otherwise, the algorithm continues with the next iteration step. Each

319

step consists of (i) solving the discretized problem (Section 3.1) and (ii) solving

320

NLPℓ for all ℓ ∈ L with parameters and bounds computed by (27) – (29) to

321

obtain a primal MINLP solution for the original problem (Section 3.2). After

322

that and in the case of nontermination, (iii) the discretization is adapted as it

323

is described in Section 3.3. All parameters concerning the discretization grid,

324

val val val including the operation discretizations V˙ skℓ1 < V˙ skℓ2 < . . . < V˙ skℓj max , which

cr

us

an

M

skℓ

are based on the discrete sizes, and the formulation of the discretized problem

Ac ce p

te

d

325

ip t

311

Figure 5: Adaptive Discretization Algorithm

18

Page 19 of 39

(16) – (26) are updated. The next iteration step of the algorithm starts with

327

solving the adapted discretized problem. Notice that the last iteration’s solution

328

of the discretized problem can be used to warm-start the discretized MIP, since

329

the discrete decisions are still feasible in the adapted grid. This improves the

330

performance of solving the MIP.

331

Remark 1. Alternative Approximations for the Discretized Problem. To en-

332

able equality in the energy balances in the discretized problem, we expand the

333

pure discretization of the operation to a piecewise linearization (Section 3.1).

334

However, this is not necessary for the functionality of the proposed algorithm,

335

since the discretized MIP is only used as an approximation of the original prob-

336

lem. We developed and tested other approaches for approximating the problem

337

via discretization in addition to the discretized MIP (16) – (26). In a pure

338

discretization, without the piecewise linearization, we can only request for at

339

least fulfilling the energy demands. Consequently, energy excesses are possi-

340

ble. These excesses can occur since it is more favorable to produce electricity

341

using CHP engines than to purchase it from the power grid. Thus, additional

342

heat production by CHP engines can decrease operational costs of DESS. For

343

that reason, there is more energy excess in solutions of a pure discretization

344

than might be expected. Consequently, the approximation quality gets worse

345

since overproduction is not allowed in the original problem. We develop two

346

approaches to restrict the amount of energy excess in the pure discretization

347

and to improve the approximation quality: First, by adding new and nontrivial

348

constraints to the discretized problem, secondly, by penalizing energy excess in

349

the objective function. It turns out that our adaptive discretization algorithm

350

using the discretized problem with piecewise linearization of Section 3.1 outper-

351

forms these more sophisticated alternatives on the set of test instances (Section

352

4.1).

Ac ce p

te

d

M

an

us

cr

ip t

326

19

Page 20 of 39

4. Computational Study and Results

354

In the computational study, we analyze the solution quality of our adaptive

355

discretization algorithm (short AdaptDiscAlgo, Section 3) in comparison to (i)

356

primal solutions of MINLP (Section 2.2) computed with BARON and to (ii)

357

approximate solutions of piecewise linearized models following the explanations

358

in Section 2.3. Furthermore, we examine the running times of all solving ap-

359

proaches.

360

Section 4.1 gives an overview and references of the online available considered

361

test instances and Section 4.2 contains some details on the implementation of

362

all approaches for computing MINLP solutions to synthesis of DESS. The com-

363

putational results are presented in Section 4.3.

364

4.1. Problem Instances

365

For the computational study on the performance our algorithm, we derived a

366

test-set: DESSLib Bahl et al. (2015). The DESSLib contains categorized prob-

367

lem instances based on the original real-world example stated by Voll et al.

368

(2013b). The categories are characterized by 2 dimensions, the number of con-

369

sidered units in the superstructure and the number of considered load cases.

370

The category for the number of considered units is denoted by S4, S8, S12, S16,

Ac ce p

te

d

M

an

us

cr

ip t

353

371

e.g., S4 corresponds to one unit for each type of equipment. The number of

372

considered load cases is denoted by L1, L2, L4, L6, L8, L12, L16, L24. Each cat-

373

egory (e.g., S8L4) consists of 10 instances with stochastic variations around

374

the original demand time-series. We assinged stochastic variations with latin

375

hypercube sampling (Iman et al., 1981) and a variation of ±5% of the original

376

demand. Thus, the resulting DESSLib consists of 320 problem instances. We

377

use the DESSLib to evaluate the performance of the proposed adaptive dis-

378

cretization algorithm. The short notation, e.g., S4L{1, 2} denotes the set of all

379

test instances of categories S4L1 and S4L2.

20

Page 21 of 39

4.2. Implementation

381

Software and Machine. Our adaptive discretization algorithm (Section 3) was

382

implemented in GAMS 24.4.3 (GAMS Development Corporation, 2015) using

383

its Python-API. Computations are performed on one core of a Linux maschine

384

with 3.40GHz Intel Core i7-3770 processor and 32 GB RAM. To solve mixed-

385

integer linear programs, i.e., discretized MIP (16) – (26) and MILP by Voll et al.

386

(2013b) (Section 2.3), we use the default setting of CPLEX 12.6.1.0 (IBM, 2015).

387

To solve (mixed-integer) nonlinear programs, i.e., MINLP (1) – (12), feasibility

388

NLP (Section 3.2) and MINLPlin,feas (13) – (15), we use the default setting

389

of BARON 14.4.0 (Tawarmalani & Sahinidis, 2005). BARON was selected as

390

MINLP solver due its robustness in a preliminary study.

391

Discretization Grid. The number of discrete sizes ksmax and the number of dis-

392

max crete operations jskℓ in the discretized problem (Section 3.1) are ksmax = 5 and,

393

max with some exceptions, jskℓ = 10. However, if the interval of possible operations

394

N,val ˙ N,val max [αsmin V˙ sk , Vsk ] is smaller than 1800 kW, we reduce jskℓ ≤ 10 as much as

395

val,diff max necessary so that V˙ skℓj ≥ 200 kW or jskℓ = 2 holds. These parameters have

396

been determined in preliminary studies.

397

Limits and Stop Criterion. For solving the discretized MIP (Section 3.1), a

Ac ce p

te

d

M

an

us

cr

ip t

380

398

time limit of 300 seconds is implemented in the very first iteration step of an

399

algorithm run. For any further step, this time limit is set to 100 seconds and the

400

last iteration’s solution is used as a starting point. The time limit for solving

401

the feasibility NLP (Section 3.2) is set to 100 seconds. Limits for the optimality

402

gap are 0.1% (discretized MIP) and 0.001% (feasibility NLP). The adaptive

403

discretization algorithm terminates, if the running time reaches the limit of one

404

hour or the improvement of the best MINLP solution is less than 0.1% over the

405

last two iteration steps.

406

For the benchmark MINLP (1) – (12), the same limit of one hour running time

407

is implemented in each case. To obtain a feasible MINLP solution based on an

408

approximate MILP solution, we solve MILP by Voll et al. (2013b) (Section 2.3)

21

Page 22 of 39

with a time limit of one hour and thereafter we solve MINLPlin,feas (13) – (15)

410

with a time limit of one hour as well.

411

Preprocessing on Discretization Grid. Depending on the input data, particu-

412

N ˙ larly the demands of the load cases, some binary variables V˙ sk , Vskℓj can be

413

eliminated by preprocessing. For chiller units s ∈ A ∪ T , the operation variables

414

val V˙ skℓj ∈ {0, 1} with V˙ skℓj > E˙ ℓcool , i.e., supply greater than demand, will not be

415

part of any feasible solution of the discretized problem. In analogy, for heat-

416

val producing units s ∈ B∪C, an upper bound for V˙ skℓj is given by the heat demand

417

E˙ ℓheat plus the maximal possible heat demand of absorption chillers s ∈ A. If

418

for indices s, k these constraints remove all discrete operation variables V˙ skℓj ,

419

N can be eliminated as well. This the corresponding discrete size variable V˙ sk

420

preprocessing on the discretization grid is quite natural, however, its effect may

421

not be underestimated, because of the large number of binary variables in the

422

discretized problem.

423

4.3. Computational Results

424

Before we analyze the robustness of our proposed algorithm compared to the two

425

benchmarking approaches on basis of the entire set of 320 instances, we focus

426

on one randomly chosen instance and compare the numerical solution results,

Ac ce p

te

d

M

an

us

cr

ip t

409

427

i.e., DESS structure and equipment sizes, of all three considered approaches.

428

Table 1 shows the DESS structure and equipment sizes (numbers rounded) for

429

DESSLib’s instance S8L4 8 obtained by the solution approaches MINLP (Sec-

430

tion 2.2), MINLPlin,feas (Section 2.3), and AdaptDiscAlgo (Section 3). Table

431

1 does not contain the DESS resulting from the MILP of Voll et al. (2013b),

432

which was the basis for the MINLPlin,feas solution (cf. Section 2.3). The MILP

433

solution differs from the MINLPlin,feas solution just in a small increase of the

434

boiler’s sizing: from 11.45 MW to 12.02 MW. By this change, the solution

435

becomes feasible. The three approaches mentioned in Table 1 lead to three

436

different DESS structures, whereby the solution of AdaptDiscAlgo has the best

437

objective value, i.e., net present value. In fact, AdaptDiscAlgo’s solution has,

22

Page 23 of 39

MINLP

MINLPlin,feas

AdaptDiscAlgo

Boiler #1

12.02 MW

11.45 MW

10.53 MW

Boiler #2

-

-

0.10 MW

CHP engine #1

3.20 MW

2.50 MW

3.20 MW

CHP engine #2

2.35 MW

2.30 MW

2.15 MW

Turbo chiller #1

7.88 MW

3.17 MW

Turbo chiller #2

-

1.88 MW

Absorption chiller #1

6.05 MW

6.43 MW

6.50 MW

-

2.43 MW

2.07 MW

cr 3.48 MW 1.90 MW

us

−6.81 · 10

Net present value

7

−5.00 · 10

7

an

Absorption chiller #2

ip t

(Instance S8L4 8)

−4.85 · 107

Table 1: DESS structure and equipment sizes computed by solution approaches MINLP,

M

MINLPlin,feas , and AdaptDiscAlgo for DESSLib’s instance S8L4 8.

compared to the computed result of MINLP, a 28% higher net present value.

439

With regard to MINLPlin,feas , the presented AdaptDiscAlgo leads to a slightly,

440

i.e., 3%, better solution. Looking at the DESS structures it turns out that using

441

a second compressor chiller and a second absorption chiller is profitable in this

442

problem instance. Furthermore, the total size of the boilers decrease, while the

te

d

438

net present value of the three DESS systems increases.

Ac ce p 443

444

445

In the further course of the computational results we evaluate the performance

446

of AdaptDiscAlgo regarding objective value and computation times on the basis

447

of the entire set of 320 problem instances. Among other things, we will see that

448

the instance chosen for Table 1 is a good representative, in the sense that the

449

solutions of AdaptDiscAlgo outperform the two benchmarking approaches.

450

451

For evaluating the quality of a feasible solution of a problem, we consider the

452

primal gap to the objective value of a given reference solution, i.e., an optimal

453

or other known solution.

23

Page 24 of 39

454

Definition 1. The primal gap of a feasible solution x with objective value f (x)

455

and a reference solution x⋆ with objective value f (x⋆ ) is, except for trivial cases,

456

defined by gapp (x, x⋆ ) :=

457

We compute averages over parts of the test instances using the shifted geometric

458

mean, which is customary in computational optimization, see Achterberg (2007).

459

Definition 2. The shifted geometric mean of numbers a1 , . . . , ak ∈ R and a

460

shift ζ ∈ R+ with (ai + ζ) > 0, i = 1, . . . , k is defined by γζ (a1 , . . . , ak ) :=

k ∏

) k1 (ai + ζ)

i=1

− ζ.

(30)

an

(

us

cr

ip t

f (x)−f (x⋆ ) |f (x⋆ )| .

We use a shift of ζ = 100 for time in seconds and values in percent, e.g., primal-

462

dual bound, and a shift of ζ = 10 for number of iterations.

M

461

463

The computational results of the solution quality of the benchmark MINLP

465

and AdaptDiscAlgo is represented by a heatmap in Figure 6. For every test in-

466

stance category, the average of relative improvement of the best solution of

467

AdaptDiscAlgo in comparison to the best solution of MINLP (1) – (12) is

468

indicated by the percentage and coloring of the heatmap square. To put it

te

d

464

briefly, the greener, the better the solution quality of AdaptDiscAlgo compared

470

to MINLP. The averages are calculated with the best MINLP solution as ref-

471

erence solution (Definition 1) and only instances are considered where both

472

approaches provide an MINLP solution. Additionally and enclosed in brack-

473

ets, the number of test instances (max. ten per category) is given, for which

474

AdaptDiscAlgo (right number) respectively MINLP (left) computes an MINLP

475

solution within the time limit.

476

For the smallest test instances, MINLP was able to solve every test instance of

477

categories S4L{1, 2, 4, 6} optimally within the time limit. For 63% (202 out of

478

320) of all considered test instances, at least a feasible solution was computed

479

by MINLP, in all other cases the time limit was reached without any primal

480

solution. In contrast, AdaptDiscAlgo was able to compute an MINLP feasible

Ac ce p 469

24

Page 25 of 39

solution for every single test instance. Of course, the proposed algorithm can

482

not ensure optimality of its solutions. However, AdaptDiscAlgo provides near-

483

optimal solutions for test instances solved optimally by MINLP. Irritatingly, in

484

category S4L6, AdaptDiscAlgo computes on an average a 0.1%-better solution

485

than an (according to BARON) globally optimal solution of MINLP. This is due

486

to the fact that BARON cuts off optimal solutions during the solution process

Ac ce p

te

d

M

an

us

cr

ip t

481

Figure 6: Heatmap on the improvement of solution quality of AdaptDiscAlgo in relation to MINLP (1) – (12); averaged for each test instance category. In brackets: number of test instances (max. 10 per category), for which MINLP (left) respectively AdaptDiscAlgo (right) computes at least one MINLP solution in the time limit.

25

Page 26 of 39

for some of these instances. Unfortunately, this is not unusual in computational

488

(nonlinear) optimization due to limited machine accuracy. All in all, except for

489

the smallest instance categories, AdaptDiscAlgo is able to compute up to 40%

490

better MINLP solutions than MINLP.

491

All primal bounds, dual bounds and computation times of the MINLP and

492

AdaptDiscAlgo of the computations executed for this work are online available

493

at DESSLib (Bahl et al., 2015).

us

cr

ip t

487

Ac ce p

te

d

M

an

494

Figure 7: Heatmap on the improvement of solution quality of AdaptDiscAlgo in relation to MINLPlin,feas ; averaged for each test instance category.

26

Page 27 of 39

Since common DESS-solving approaches consist of piecewise linearized models

496

(Section 1), we compare approximate solutions of such a model with MINLP

497

solutions of the proposed AdaptDiscAlgo following the explanations in Section

498

2.3. The structure ys∗ and sizing decisions VsN∗ of MILP solutions are in 85%,

499

i.e., 273 out of 320 test instances, not feasible in the original MINLP problem. If

500

it was feasible, it was a test instance of small-sized category S4. For this reason,

501

we consult MINLPlin,feas (13) – (15) to compute a MINLP-feasible solution close

502

to the approximate MILP solution and evaluate it with the original objective (1).

503

The comparison of the solution quality of MINLPlin,feas and AdaptDiscAlgo is

504

shown by a heatmap in Figure 7. The averages are calculated with MINLPlin,feas

505

solutions as reference solutions (Definition 1). In all categories, AdaptDiscAlgo

506

provides comparably good or slightly better solutions as post-processed solutions

507

of the piecewise-linearized approach.

Ac ce p

te

d

M

an

us

cr

ip t

495

Figure 8: Average running times of MINLP, MILP, MINLPlin,feas and AdaptDiscAlgo.

508

The running times of all approaches are shown in Figure 8. Since only instances

509

of categories S4L{1, 2, 4, 6} are solved optimally by MINLP, for all other test in-

510

stances, the running time reaches the limit of one hour. The piecewise linearized

511

model (MILP) runs for most instances less than one hour. However, to obtain

512

an MINLP-feasible solution through MINLPlin,feas , this whole approach runs up

513

to two hours. In comparison, the proposed AdaptDiscAlgo terminates generally

514

after a fraction of an hour and thus outperforms the compared approaches in 27

Page 28 of 39

terms of solution quality as well as running time, except for the very smallest

516

test instances.

517

The running time of AdaptDiscAlgo is split up over the parts of the algorithm as

518

follows: On average 93% of the running time is used for solving the discretized

519

MIP (Section 3.1) and 5 iteration steps are passed on an average until the stop

520

criterion of convergence is satisfied. For all considered test instances, the fea-

521

sibility NLP (Section 3.2) was feasible in every single iteration step. This was

522

not the case for the alternative approaches of the discretized problem which are

523

briefly mentioned in Section 3.4.

us

cr

ip t

515

an

524

5. Summary and Outlook

526

The superstructure-based synthesis of decentralized energy supply systems can

527

be formulated as a mixed-integer nonlinear program. By including, e.g., nonlin-

528

ear part-load performances, investment costs, and strict energy balances, this

529

optimization problem is unavoidably nonconvex. In this paper, we do not cir-

530

cumvent this issue via approximating the problem by linearizing the perfor-

531

mance and investment cost models of all considered types of equipment. Such

532

approximate solutions from linearized problems are not necessarily feasible for

533

the original nonlinear problem. In this paper, we set the focus on computing

534

solutions which are feasible for the nonlinear synthesis problem. Since global

535

MINLP solvers, e.g., BARON, have difficulties to provide primal solutions for

536

real-world-based test instances, we propose a problem-specific solution approach

537

for the nonlinear synthesis of DESS. This optimization-based algorithm consists

538

of a discretized and linearly formulated version of the synthesis problem, whose

539

underlying discretization grid is iteratively adapted. The resulting approximate

540

problem is of such a nature that solution-specific decisions are also feasible in

541

the original nonlinear problem and by solving a decomposable NLP, this leads

542

to MINLP solutions. A computational study based on a set of test instances

543

obtained from real industrial data shows that the proposed adaptive discretiza-

Ac ce p

te

d

M

525

28

Page 29 of 39

tion algorithm computes better MINLP solutions in less computation times than

545

state-of-the-art solvers. Thus, the proposed algorithm provides an efficient so-

546

lution method to the synthesis of decentralized energy supply systems.

ip t

544

547

In future work, one should expand the adaptive discretization algorithm con-

549

cerning methods to compute dual bounds for the MINLP to estimate the op-

550

timality gap of the algorithm’s primal solutions. The (mostly very weak) dual

551

bounds provided by BARON does not contribute meaningful information.

552

A mathematically feasible or even optimal solution is usually only an approxi-

553

mation of a real-world implementation, since a model never represents the real

554

problem perfectly. Real world decisions might be influenced by constraints not

555

represented in the model, e.g., (missing) maintenance knowledge in the com-

556

pany for specific technologies. In Voll et al. (2013a) and Hennen et al. (2016)

557

we show that several near-optimal solution alternatives exist. Analyzing these

558

near-optimal solutions allows to derive real-world decision options. Since the

559

proposed AdaptDiscAlgo efficiently provides feasible solutions of the nonlinear

560

synthesis problem, in future work the algorithm could be expanded to efficiently

561

generate many near-optimal solution alternatives.

562

To identify the suitable superstructure, one could complement the proposed

Ac ce p

te

d

M

an

us

cr

548

563

algorithm with the successive superstructure expansion method of Voll et al.

564

(2013b).

565

The application of the proposed discretization algorithm to other hard to solve

566

synthesis problems seems quite promising.

29

Page 30 of 39

References

568

Achterberg, T. (2007). Constraint Integer Programming. Ph.D. thesis, Techni-

569

ip t

567

sche Universit¨ at Berlin.

Arcuri, P., Florio, G., & Fragiacomo, P. (2007). A mixed integer programming

571

model for optimal design of trigeneration in a hospital complex. Energy, 32 ,

572

1430–1447. doi:10.1016/j.energy.2006.10.023.

574

us

573

cr

570

ASUE (2011). BHKW-Kenndaten 2011. In Datenblatt. Arbeitsgemeinschaft f¨ ur sparsamen und umweltfreundlichen Energieverbrauch e.V.

Bahl, B., Goderbauer, S., Arnold, F., Voll, P., Bardow, A., Koster, A., &

576

L¨ ubbecke, M. (2015). DESSLib — Test Instances for Optimization of Decen-

577

tralized Energy Supply Systems. URL: http://www.math2.rwth-aachen.

578

de/DESSLib [Online; accessed 07-June-2016].

M

579

an

575

Bruno, J., Fernandez, F., Castells, F., & Grossmann, I. (1998).

A

rigorous MINLP model for the optimal synthesis and operation of

581

utility plants.

582

//www.sciencedirect.com/science/article/pii/S0263876298716430.

583

doi:10.1205/026387698524901.

d

580

Ac ce p

te

Chem. Eng. Res. Des., 76 , 246–258. URL: http:

584

Chen, C.-L., & Lin, C.-Y. (2011).

A flexible structural and operational

585

design of steam systems. Appl. Therm. Eng., 31 , 2084–2093. URL: http:

586

//www.sciencedirect.com/science/article/pii/S1359431111001396.

587

doi:10.1016/j.applthermaleng.2011.03.007.

588

Dimopoulos, G. G., & Frangopoulos, C. A. (2008).

Optimization of en-

589

ergy systems based on evolutionary and social metaphors.

590

33 , 171–179. URL: http://www.sciencedirect.com/science/article/

591

pii/S0360544207001600. doi:10.1016/j.energy.2007.09.002.

592

Drumm,

C.,

Busch,

J.,

Dietrich,

W.,

Eickmans,

J.,

Energy,

& Jupke,

A.

STRUCTese ® - Energy efficiency management for the pro-

593

(2013).

594

cess industry.

Chem. Eng. Process:

Process Intensification,

67 ,

30

Page 31 of 39

595

99–110.

596

S0255270112001845. doi:10.1016/j.cep.2012.09.009.

http://www.sciencedirect.com/science/article/pii/

ip t

URL:

Fabrizio, E. (2008). Modelling of multi-energy systems in buildings. Ph.D. thesis

598

, Politecnico di Turino and Institut National des Sciences Appliqu´ees de Lyon.

599

Farkas, T., Rev, E., & Lelkes, Z. (2005). Process flowsheet superstructures:

600

Structural multiplicity and redundancy. Computers & Chemical Engineering,

601

29 , 2198–2214. URL: http://dx.doi.org/10.1016/j.compchemeng.2005.

602

07.008. doi:10.1016/j.compchemeng.2005.07.008.

us

cr

597

Frangopoulos, C. A., Spakovsky, M. R. v., & Sciubba, E. (2002). A brief re-

604

view of methods for the design and synthesis optimization of energy systems.

605

International Journal of Applied Thermodynamics, 5 , 151–160.

607

M

606

an

603

GAMS Development Corporation (2015). General Algebraic Modeling System (GAMS). URL: http://www.gams.com.

Geißler, B., Martin, A., Morsi, A., & Schewe, L. (2011). Using piecewise linear

609

functions for solving MINLPs. In Mixed Integer Nonlinear Programming (pp.

610

287–314). Springer Science + Business Media. URL: http://dx.doi.org/

611

10.1007/978-1-4614-1927-3_10. doi:10.1007/978-1-4614-1927-3_10.

Ac ce p

te

d

608

612

Gupte, A., Ahmed, S., Cheon, M. S., & Dey, S. (2013). Solving mixed integer

613

bilinear problems using MILP formulations. SIAM J. Optim., 23 , 721–744.

614

URL: http://dx.doi.org/10.1137/110836183. doi:10.1137/110836183.

615

Hennen, M., Lampe, M., Voll, P., & Bardow, A. (2016). How to explore and

616

analyze the decision space in the synthesis of energy supply systems. In Pro-

617

ceedings of ECOS 2016 - The 29th International Conference on Efficiency,

618

Cost, Optimization, Simulation and Environmental Impact of Energy Sys-

619

tems.

620

621

IBM (2015). IBM ILOG CPLEX Optimization Studio. URL: http://www. ilog.com/products/cplex [Online; accessed 7-December-2015].

31

Page 32 of 39

Iman, R. L., Campbell, J., & Helton, J. (1981). An approach to sensitivity

623

analysis of computer models. i- introduction, input, variable selection and

624

preliminary variable assessment. Journal of quality technology, 13 , 174–183.

625

IUTA (2002). Preisatlas - Ableitung von Kostenfunktionen f¨ ur Komponenten

626

der rationellen Energienutzung. In Abschlussbericht. Institut f¨ ur Energie- und

627

Umwelttechnik e.V. M.,

Fisk,

D.,

& Shah,

N. (2014).

cr

Jennings,

Modelling and op-

us

628

ip t

622

629

timization of retrofitting residential energy systems at the urban

630

scale.

631

com/science/article/pii/S0360544213009432. doi:http://dx.doi.org/

632

10.1016/j.energy.2013.10.076.

64 ,

220–233.

URL:

http://www.sciencedirect.

an

Energy,

Kolodziej, S., Castro, P. M., & Grossmann, I. E. (2013). Global optimization of

634

bilinear programs with a multiparametric disaggregation technique. Journal

635

of Global Optimization, 57 , 1039–1063. URL: http://dx.doi.org/10.1007/

636

s10898-012-0022-1. doi:10.1007/s10898-012-0022-1.

d

M

633

Leyffer, S., Sartenaer, A., & Wanufelle, E. (2008). Branch-and-refine for mixed-

638

integer nonconvex global optimization. Preprint ANL/MCS-P1547-0908, Ar-

te

637

gonne National Laboratory, Mathematics and Computer Science Division.

640

Lozano, M. A., Ramos, J. C., Carvalho, M., & Serra, L. M. (2009). Structure op-

641

timization of energy supply systems in tertiary sector buildings. Energ Build ,

642

41 , 1063–1075. URL: http://www.sciencedirect.com/science/article/

643

pii/S037877880900111X. doi:10.1016/j.enbuild.2009.05.008.

Ac ce p 639

644

Mar´echal, F., & Kalitventzeff, B. (2003). Targeting the integration of multi-

645

period utility systems for site scale process integration. Appl. Therm. Eng.,

646

23 , 1763–1784. URL: http://www.sciencedirect.com/science/article/

647

pii/S135943110300142X. doi:10.1016/S1359-4311(03)00142-X.

648

Mar´echal, F., Weber, C., & Favrat, D. (2008). Multiobjective design and opti-

649

mization of urban energy systems. In E. N. Pistikopoulos, M. C. Georgiadis, 32

Page 33 of 39

E. S. Kikkinides, & V. Dua (Eds.), Energy systems engineering. Weinheim,

651

Germany: Wiley-VCH volume 5 of Process systems engineering.

652

Østergaard,

P. A. (2009).

ip t

650

Reviewing optimisation criteria for en-

ergy systems analyses of renewable energy integration.

Energy, 34 ,

654

1236–1245. URL: http://www.sciencedirect.com/science/article/pii/

655

S0360544209001777. doi:10.1016/j.energy.2009.05.004.

cr

653

Papalexandri, K. P., Pistikopoulos, E. N., & Kalitventzeff, B. (1998). Modelling

657

and optimization aspects in energy management and plant operation with

658

variable energy demands-application to industrial problems. Comput. Chem.

659

Eng., 22 , 1319–1333. URL: http://www.sciencedirect.com/science/

660

article/B6TFT-3WTNY92-J/2/c4619100d5e2c4581678da2f2b88e278.

661

doi:DOI:10.1016/S0098-1354(98)00016-7.

662

Papoulias, S. A., & Grossmann, I. E. (1983).

M

an

us

656

A structural optimization

663

approach in process synthesis - I: Utility systems.

664

Eng.,

665

article/B6TFT-43PRJJC-20/2/d103399965d19f265da6d82a661cb8a9.

666

doi:DOI:10.1016/0098-1354(83)85022-4.

http://www.sciencedirect.com/science/

d

URL:

te

695–706.

Pham, V., Laird, C., & El-Halwagi, M. (2009). Convex hull discretization

Ac ce p

667

7,

Comput. Chem.

668

approach to the global optimization of pooling problems. Industrial & Engi-

669

neering Chemistry Research, 48 , 1973–1979. URL: http://dx.doi.org/10.

670

1021/ie8003573. doi:10.1021/ie8003573.

671

Prokopakis, G. J., & Maroulis, Z. B. (1996). Real-time management and opti-

672

misation of industrial utilities systems. Comput. Chem. Eng., 20, Supplement

673

1 , S623 – S628. URL: http://www.sciencedirect.com/science/article/

674

pii/0098135496001135. doi:10.1016/0098-1354(96)00113-5.

675

Stojiljkovic, M., Stojiljkovic, M., & Blagojevic, B. (2014).

Multi-objective

676

combinatorial optimization of trigeneration plants based on metaheuristics.

677

Energies, 7 , 8554–8581. URL: http://dx.doi.org/10.3390/en7128554.

678

doi:10.3390/en7128554. 33

Page 34 of 39

679

Tawarmalani, M., & Sahinidis, N. V. (2005).

A polyhedral branch-and-

cut approach to global optimization.

Mathematical Programming, 103 ,

681

225–249. URL: http://dx.doi.org/10.1007/s10107-005-0581-8. doi:10.

682

1007/s10107-005-0581-8.

ip t

680

Tong, C., Palazoglu, A., El-Farra, N. H., & Yan, X. (2015). Energy demand

684

management for process systems through production scheduling and control.

685

AIChE Journal , 61 , 3756–3769. URL: http://dx.doi.org/10.1002/aic.

686

15033. doi:10.1002/aic.15033.

us

Varbanov, P., Doyle, S., & Smith, R. (2004). Modelling and optimization

an

687

cr

683

of utility systems.

Chem. Eng. Res. Des., 82 , 561–578. URL: http:

689

//www.sciencedirect.com/science/article/pii/S0263876204725289.

690

doi:10.1205/026387604323142603.

M

688

Varbanov, P. S., Perry, S., Klemes, J. J., & Smith, R. (2005). Synthesis of indus-

692

trial utility systems: cost-effective de-carbonisation. Appl. Therm. Eng., 25 ,

693

985–1001. URL: http://www.sciencedirect.com/science/article/pii/

694

S1359431104002303. doi:10.1016/j.applthermaleng.2004.06.023.

te

d

691

Voll, P., Hennen, M., Klaffke, C., Lampe, M., & Bardow, A. (2013a). Exploring

696

the near-optimal solution space for the synthesis of distributed energy supply

Ac ce p

695

697

systems. Chemical Engineering Transactions, 35 , 277–282. URL: http://

698

doi.org/10.3303/CET1335046. doi:10.3303/CET1335046.

699

Voll, P., Klaffke, C., Hennen, M., & Bardow, A. (2013b).

Automated

700

superstructure-based synthesis and optimization of distributed energy sup-

701

ply systems. Energy, 50 , 374–388. URL: http://www.sciencedirect.com/

702

science/article/pii/S0360544212008079. doi:10.1016/j.energy.2012.

703

10.045.

704

Voll, P., Lampe, M., Wrobel, G., & Bardow, A. (2012). Superstructure-free

705

synthesis and optimization of distributed industrial energy supply systems.

706

Energy, 45 , 424–435. URL: http://www.sciencedirect.com/science/

707

article/pii/S0360544212000461. doi:10.1016/j.energy.2012.01.041. 34

Page 35 of 39

Yokoyama, R., Shinano, Y., Taniguchi, S., Ohkura, M., & Wakui, T. (2015).

709

Optimization of energy supply systems by MILP branch and bound method

710

in consideration of hierarchical relationship between design and operation.

711

Energy Convers Manage, 92 , 92 – 104. URL: http://www.sciencedirect.

712

com/science/article/pii/S0196890414010498. doi:http://dx.doi.org/

713

10.1016/j.enconman.2014.12.020.

cr

Yue, D., & You, F. (2014).

Game-theoretic modeling and optimization of

us

714

ip t

708

multi-echelon supply chain design and operation under stackelberg game

716

and market equilibrium. Computers & Chemical Engineering, 71 , 347–361.

717

URL: http://dx.doi.org/10.1016/j.compchemeng.2014.08.010. doi:10.

718

1016/j.compchemeng.2014.08.010.

an

715

Zhao, H., Rong, G., & Feng, Y. (2015). Effective solution approach for integrated

720

optimization models of refinery production and utility system. Industrial &

721

Engineering Chemistry Research, 54 , 9238–9250. URL: http://dx.doi.org/

722

10.1021/acs.iecr.5b00713. doi:10.1021/acs.iecr.5b00713.

Ac ce p

te

d

M

719

35

Page 36 of 39

Appendix A. Economic Parameters and Equipment Models

724

The economic parameters for the objective function, i.e., net present value, are

725

taken from Voll et al. (2013b) and listed in table A.2.

726

We list parameters of the considered types of equipment in A.3. Furthermore, we

727

state the nonlinear models for part-load operation and investment cost curves

728

used in this paper below. The part-load performance of CHP units is based

729

on measured data-points for several existing units. Moreover we assume that

730

the part-load operation is not depending on the size of equipment, thus scaling

731

to a normalized output power is possible. The part-load efficiency for boilers

732

and absorption chillers is modeled in analogy to Fabrizio (2008). The part-load

733

performance behavior is modeled in analogy to Fabrizio (2008) and additional

734

correspondence with turbo compression manufacturers. The nominal efficiency

735

of the CHP engines was taken from ASUE (2011). Maintenance-cost is based on

736

IUTA (2002), the investment cost curves consider are composed on information

737

from IUTA (2002) and databases of industrial partners. pel,sell

pgas,buy

i

γ CF

0.16 ct/kWh

0.10 ct/kWh

0.06 ct/kWh

0.08

10 a

Ac ce p

te

pel,buy

d

M

an

us

cr

ip t

723

Table A.2: Economic parameters of DESS synthesis problem

V˙ sN,min

V˙ sN,max

ms

αsmin

Boiler s ∈ B

0.1 MW

14 MW

1.5

0.2

CHP engine s ∈ C

0.5 MW

3.2 MW

10

0.5

Absorption chiller s ∈ A

0.05 MW

6.5 MW

1

0.2

Turbo chiller s ∈ T

0.4 MW

10 MW

4

0.2

Table A.3: Size ranges, maintenance cost factors, and minimum part-load factors of considered types of equipment.

36

Page 37 of 39

Part-load performance: (A.1) – (A.5) s ∈ B (Boiler)

(

1

U˙ s (V˙ sℓ , V˙ sN ) =

η N,B

V˙ 2 C1B · sℓ + C2B · V˙ sℓ + C3B · V˙ sN V˙ sN

)

1 = COPN,A

U˙ s (V˙ sℓ , V˙ sN )

( C1A

)

us

739

s ∈ A (Absorption chiller)

(A.1)

cr

η N,B = 0.9, C1B = 0.1021, C2B = 0.8355, C3B = 0.0666

ip t

738

V˙ 2 · sℓ + C2A · V˙ sℓ + C3A · V˙ sN V˙ sN

(A.2)

an

COPN,A = 0.67, C1A = 0.8333, C2A = −0.0833, C3A = 0.25

s ∈ T (Turbo chiller) U˙ s (V˙ sℓ , V˙ sN )

1 = COPN,T

(

M

740

C1T

V˙ 2 · sℓ + C2T · V˙ sℓ + C3T · V˙ sN V˙ sN

) (A.3)

d

COPN,T = 5.54, C1T = 0.8119, C2T = −0.1688, C3T = 0.3392

te

741

s ∈ C (CHP engine)

Ac ce p

U˙ s (V˙ sℓ , V˙ sN ) = C1C

+

C2C

V˙ sℓ · + C3C · V˙ sN + C4C · V˙ sN

(

V˙ sℓ V˙ sN

)2

( )2 + C5C · V˙ sℓ + C6C · V˙ sN

(A.4)

C1C = 550.3, C2C = −1328, C3C = −0.4537,

C4C = 668.3, C5C = 2.649, C6C = 9.571e − 05

V˙ sel (V˙ sℓ , V˙ sN ) =

V˙ sℓ C C7C + C8C · + C9C · V˙ sN + C10 · V˙ sN

(

V˙ sℓ V˙ sN

)2

( )2 C C + C11 · V˙ sℓ + C12 · V˙ sN

(A.5)

C7C = 518.8, C8C = −1203, C9C = −0.5361, C C C C10 = 579.3, C11 = 1.464, C12 = 7.728e − 05 742

37

Page 38 of 39

743

Investment cost: (A.6) – (A.9)

ip t

s ∈ B (Boiler)

cr

I(V˙ sN ) = ( [( ) )] 1.85484· 11418.6 + 64.115 · V˙ sN 0.7978 · 1.046 · 1.0917 − 1.1921 · 10−6 · V˙ sN

(A.6)

us

744

s ∈ A (Absorption chiller)

745

s ∈ T (Turbo chiller)

M

( ) I(V˙ sN ) = 0.8102 · V˙ sN · 179.63 + 4991.3436 · V˙ sN −0.6794

te

d

746

s ∈ C (CHP engine) I(V˙ sN )

(A.7)

an

I(V˙ sN ) = 0.50401 · 17554, 18 · V˙ sN 0.4345

= 9332.6 ·

Ac ce p

ηsN,th (V˙ sN ) = 0.498 − 3.55 · 10−5 · V˙ sN ,

(

V˙ sN

ηsN,el (V˙ sN ) · N,th ηs (V˙ sN )

(A.8)

)0.539

ηsN,el (V˙ sN ) = ηsN − ηsN,th (V˙ sN ),

(A.9) ηsN = 0.87

38

Page 39 of 39