arXiv:1412.4504v1 [math.

OC] 15 Dec 2014

Implementing a Unit Commitment Power Market

Model in FICO Xpress-Mosel1

René Brandenberg Matthias Silbernagl

Technische Universität München Technische Universität München

1
Preprint; Published version available at http://community.fico.com/docs/DOC-1365
 Abstract

Starting from a basic Unit Commitment problem as published in [3], we develop a fully functional,
running implementation for Xpress Mosel, useable in power market modeling. The constraints of
the model are discussed in details and solutions to left open implementation problems are presented,
guidelines on how to handle data input (from both, databases and Excel), as well as data output (to
Excel) are given.
CONTENTS I

Contents
1 Introduction 1

2 The Unit Commitment Problem 2

2.1 Basic notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.2 Physical Constraints of the Units . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.2.1 Unit Variables (constraints 2.2 and 2.3 in subsection 2.7) . . . . . . . . . 2
2.2.2 Minimal Up- and Downtime (constraints 2.7 to 2.10 in subsection 2.7) . . 3
2.2.3 Minimal and Maximal Production (constraint 2.11 in subsection 2.7) . . . 3
2.2.4 Ramping (constraints 2.12 to 2.14 in subsection 2.7) . . . . . . . . . . . 4
2.2.5 Storage (constraints 2.4, 2.15 to 2.19 in subsection 2.7) . . . . . . . . . . 5
2.3 Power Grid Constraints (constraints 2.20 and 2.21 in subsection 2.7) . . . . . . . 6
2.4 Cost Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.4.1 Objective Function (constraints 2.1 and 2.5 in subsection 2.7) . . . . . . 6
2.4.2 Production Costs (constraint 2.22 in subsection 2.7) . . . . . . . . . . . . 6
2.4.3 Startup and Shutdown Costs (constraints 2.23 and 2.24 in subsection 2.7) 7
2.4.4 Thinning Out the Startup Cost Function . . . . . . . . . . . . . . . . . . 8
2.5 Dealing with Infeasibilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.6 Parameters summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.6.1 General parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.6.2 Parameters for each period . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.6.3 Parameters for each power unit . . . . . . . . . . . . . . . . . . . . . . . 11
2.7 Model summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3 Implementing the Model 14

3.1 Parameters (subsection 2.6) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.2 Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.3 Physical Constraints of the Units (constraints 2.7 to 2.14) . . . . . . . . . . . . . 16
3.4 Power Grid Constraints (constraints 2.20 and 2.21) . . . . . . . . . . . . . . . . 17
3.5 Production Cost Constraints (constraint 2.22) . . . . . . . . . . . . . . . . . . . 17
3.6 Startup and Shutdown Costs (constraints 2.23 and 2.24) . . . . . . . . . . . . . . 17
3.6.1 Thinning Out the Startup Cost Function . . . . . . . . . . . . . . . . . . 17
3.7 Objective Function (constraint 2.1) . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.8 Dealing with Infeasibilities (subsection 2.5) . . . . . . . . . . . . . . . . . . . . 19

4 Data Input 20
4.1 General Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.2 Unit Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.3 Period Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
CONTENTS II

5 Data Output 25
5.1 Electricity price . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
5.2 Postprocessing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
5.3 Data Output to Excel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

6 Remarks 29

A Excel Functions 31
A.1 Function isExcelDocument . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
A.2 Referencing cells in Excel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
A.2.1 Column index conversions in the A1 notation . . . . . . . . . . . . . . . 32
A.3 Function Join . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
1 INTRODUCTION 1

1 Introduction
The deregulation of the energy market in the last two decades attracted the researchers’ attention,
as electricity producers constantly needed to optimize their production to stay competitive. The
raising volatility in energy production, mostly due to the renewable resources, even increases the
potential for optimization and therefore the interest in power market models.
We consider a Unit Commitment (UC) model based on the mixed integer linear program
(MILP) presented in [3], which assumes perfect competition, and as such approximates the typical
oligopolistic power market in deregulated countries. Its main advantage over possibly more
sophisticated game theoretical or stochastic approaches is its capability to handle real-world
problem sizes with an accurate physical model of the power units, combined with its relatively
good comprehensibility.
To make the model in [3] suitable for real-life application, we need to modify it in a few key
aspects. In particular, we focus on

• the reduction of the prohibitively large number of constraints needed to model the start-up
costs,

• the incorporation of storage units,

• the analysis of causes of infeasibility and how to avoid them, and

• the incorporation of daily fuel prices, as fuel costs make up the better part of the production
costs.

After a short theoretical presentation of the model, we put our main focus on a discussion of
the constraints and a thorough guide on their implementation.
Finally, we show how to read the needed input data (which may be considerably huge, if long
time period real world data is used) from either an SQL database or an Excel spreadsheet. The
output of the results is written to an Excel spreadsheet, for further processing and analysis.
We use the Xpress-Mosel modeling language [1], which is part of the Xpress Optimization
Suite [2].
2 THE UNIT COMMITMENT PROBLEM 2

2 The Unit Commitment Problem

In the Unit Commitment problem one has to satisfy a certain energy demand over a number of
periods, using a set of given power units, while minimizing the total incurred costs.
A typical Unit Commitment problem can roughly be split into three parts:

1. The modeling of the physical constraints of each unit (subsection 2.2)

In this part the technical limitations of each individual unit are considered. We include

• the feasible range for the production level,

• the feasible range for the change rate of the production level,
• the minimal up- and downtime and
• the storage constraints.

2. The modeling of the power grid (subsection 2.3)

A model of the undelying power grid may be considered, too, which we exclude in this paper.
Hence, the overall production is just connected to the demand, a typical assumption in many
models, usually described by the keyword “big copper plate”.

3. The modeling of the costs (subsection 2.4)

Besides the obvious fuel costs, we include start up and shutdown costs and some penalties.

In the remaining of this section we discuss the components of our model in the above listed order
(subsections 2.2 - 2.4), give a short note on how to handle infeasibilities 2.5, and include summaries
of the parameters 2.6 and the model 2.7. The latter may be used as a convenient reference in the
following sections, since it is cross-referenced with the discussion and implementation over there.

2.1 Basic notation

We denote the set of time periods as K = {1, . . . , T }, and the set of units as J. The remaining
parameters will be introduced “as we go”, and summarized in subsection 2.6. We use the shortcut
notations [a..b] to denote the discrete set [a, b] ∩ Z, and [b] to denote [1..b].

2.2 Physical Constraints of the Units

2.2.1 Unit Variables (constraints 2.2 and 2.3 in subsection 2.7)

In any period k ∈ K, the units j ∈ J are modeled by their operational state (on/off) vjk , their
current production pkj and their maximal possible production pkj

∀ j ∈ J, k ∈ K : vjk ∈ {0, 1}
pkj , pkj ∈ R+ ,

where vjk = 1 iff unit j is operational in period k and may produce power.
2 THE UNIT COMMITMENT PROBLEM 3

The maximal possible production pkj is used to measure the spinning reserve: In case of a power
outage, the power grid must be stabilized by ramping up the currently operational units. The
available additional production capacity pkj − pkj is called spinning reserve of unit j and a general
grid constraint regulates the minimal need of spinning reserve over all units (see subsection 2.3
below).

2.2.2 Minimal Up- and Downtime (constraints 2.7 to 2.10 in subsection 2.7)

Units cannot be started up and shut down arbitrarily. For example, for a coal unit, after being shut
down, an appropriate cool-off time has to be kept.
To model this, minimal up- and downtime parameters UTj and DTj for each unit are given,
denoting the minimal time a unit has to stay operational after a start up, and the minimal cool-off
time after a shutdown, respectively. Since we do not know when a unit has last been started up or
shut down prior to the modeled time range, we further expect to be given the two parameters IUTj
and IDTj , which tell us for how long a unit needs to be operational or shut down initially. Given
these parameters, we model the minimal uptime as

∀ j ∈ J, k ∈ [IUTj ] : vjk = 1

for the first periods, and as

∀ j ∈ J, k ∈ [(IUTj + 2)..T ], i ∈ [UTj − 1] ∩ [T − k] : vjk+i ≥ vjk − vjk−1

for the remaining periods. Consider the right-hand side vjk − vjk−1 of the latter constraint: it always
lies in {−1, 0, 1}, and is equal to 1 if and only if unit j starts up in period k. Thus, in the start-up
case, the variables vjk+i on the left-hand side are forced to 1, while the constraint is always fulfilled
otherwise. Similarly, the minimal downtime is modeled as

∀ j ∈ J, k ∈ [IDTj ] : vjk = 0

∀ j ∈ J, k ∈ [(IDTj +2)..T ], i ∈ [DTj − 1] ∩ [T − k] : vjk+i ≤ 1 − (vjk − vjk−1 )

Note that we expect IDTj , IUTj ∈ [0..T ] and DTj , UTj ∈ [T ].

2.2.3 Minimal and Maximal Production (constraint 2.11 in subsection 2.7)

If a unit j is not operational, it may not produce; otherwise both the actual and maximal possible
production level have to lie in a specific interval, defined by the two parameters Pj (minimal
production) and Pj (maximal production). All this can be expressed by

∀ j ∈ J, k ∈ K : Pj vjk ≤ pkj ≤ pkj ≤ Pj vjk

Note that storage units may have a negative minimal production Pj . However, pkj ≥ 0 only captures
power production, while power consumption is modeled by ckj (see section 2.2.5 below).
2 THE UNIT COMMITMENT PROBLEM 4

2.2.4 Ramping (constraints 2.12 to 2.14 in subsection 2.7)

Of course the production level can not change arbitrarily either. The parameters RUj / RDj denote
the maximal speed when increasing / decreasing the production level of an operating unit. For
example, a unit with RUj = 50MW/h would be able to increase it’s production level from 120MW
to 170MW in one hour, but not to 171MW. This may simply be modeled as

∀ j ∈ J, k ∈ [2..T ] : pkj ≤ pk−1

j + L · RUj
pkj ≥ pk−1
j − L · RDj vjk

where L is a parameter denoting the period length. However, at start up and shutdown units are
typically able to change their production levels faster. This higher ramping speed is denoted by the
two parameters SUj (maximal production level at start up) and SDj (maximal production level
before shutdown).
Thus for up ramping we get the constraint

∀ j ∈ J, k ∈ [2..T ] : pkj ≤ pk−1

j + L · RUj vjk−1 + SUj (1 − vjk−1 )
− min{SUj , max{Pj , 0} + L · RUj } · (1 − vjk )

The first three terms on the right-hand side ascertain that the production may only increase by
L · RUj if the unit is already running, or by SUj if the unit is starting.
The last term does not change the right hand side in the case vjk = 1, but tightens the constraint
in the case vjk = 0. Consider the constraint without the last term:

Case vjk = 0, vjk−1 = 0 : pkj ≤ pk−1

j + L · RUj vjk−1 +SUj (1 − vjk−1 ) = SUj
|{z} | {z } | {z } | {z }
=0 =0 =0 =1

Case vjk = 0, vjk−1 = 1 : pkj ≤ pk−1

j + L · RUj vjk−1 + SUj (1 − vjk−1 ) = pk−1
j + L · RUj
|{z} | {z } | {z } | {z }
=0 =1 =0

Thus, in case of vjk = 0 we could tighten the constraint by subtracting the term min{SUj , pk−1
j +
L · RUj } on the right hand side. However, since we need to multiply the term by (1 − vjk ) to leave
the case vjk = 1 unchanged, to avoid non-linearities, it should not contain a variable. Hence we
have to replace pk−1
j by its best lower bound max{Pj , 0}. This leads exactly to the last term in the
rampup constraint. This term is not necessary for a correct model, but tightens the linear relaxation,
possibly improving the solution time.
The rampdown is modeled analogously:

∀ j ∈ J, k ∈ [2..T ] : pkj ≥ pk−1

j − L · RDj vjk − SDj (1 − vjk )
+ min{SDj , max{Pj , 0} + L · RDj } · (1 − vjk−1 )
2 THE UNIT COMMITMENT PROBLEM 5

Finally, we need to put a limit on the possible maximal production in case of an imminent shutdown.
A shutdown scheduled for period k + 1 can not be postponed, even in case of a power outage.
Thus, in case of a shutdown in k + 1, the unit cannot produce more than SDj in period k:

∀ j ∈ J, k ∈ [T − 1] : pkj ≤ SDj (vjk − vjk+1 ) + Pj vjk+1

The term SDj (vjk − vjk+1 ) expresses the desired limit, while the term Pj vjk+1 keeps the constraint
valid in each of the other three settings of vjk and vjk+1 .

2.2.5 Storage (constraints 2.4, 2.15 to 2.19 in subsection 2.7)

Storage units are used to even out the demand. The classical examples for storage units are water
pumping stations and, to a lesser extent, batteries. We indicate a storage unit by a negative minimal
production Pj , meaning that the maximal storage inflow is −Pj .
Not being thermal units, the nowadays relevant storage units have very high ramping speeds, which
allow us to neglect them, i.e. we assume RUj = RDj = SDj = SUj = Pj + (−Pj ) for each
storage unit j.
In our model storage units are characterized by five parameters:

• their storage capacity SCj ,

• their storage efficiency SEj ,
• their constant inflow SIFj ,
• and their initial and final storage fill, SIj and SFj ,

and two new variables:

• the current storage fill ∀ j ∈ J, k ∈ K : skj ∈ R+ and

• the power consumption ∀ j ∈ J, k ∈ K : ckj ∈ R+ .

The maximal storage inflow and the storage capacity are now easily formulated as

∀ j ∈ J, k ∈ K : skj ≤ SCj
∀ j ∈ J, k ∈ K : ckj ≤ max{0, −Pj }

The storage fill change from period k − 1 to k has to account for the prior storage, the stored and
consumed energy, and a possible constant inflow (for example a stream feeding a reservoir):

∀ j ∈ J, k ∈ [2..T ] : skj = sk−1

j + L · (SEj ck−1
j − pk−1
j + SIFj )
2 THE UNIT COMMITMENT PROBLEM 6

The initial and final storage fill parameters SIj and SFj are set outside the model (possibly for
raising the total storage when demands and prizes are expected to increase after the considered
total time period). They apply to the storage fill at k = 1 and k = T + 1:

∀j ∈J : s1j = SIj
∀j ∈J : SFj = sTj + L · (SEj cTj − pTj + SIFj )

2.3 Power Grid Constraints (constraints 2.20 and 2.21 in subsection 2.7)
All units together have to satisfy the energy demand in each period, given as Dk :
X
∀k∈K : (pkj − ckj ) = Dk
j∈J

Moreover, to be failure-tolerant, the units should be able to compensate a power outage, caused
for example by a failing unit. This is expressed by the need to keep a given reserve Rk capacity
available:
X
∀k∈K : (pkj − pkj + ckj ) ≥ Rk
j∈J

Typically Rk is constant over time. However, a dependency e.g. on Dk maybe reasonable, too.

2.4 Cost Constraints

2.4.1 Objective Function (constraints 2.1 and 2.5 in subsection 2.7)

The objective is to minimize the overall costs consisting of production costs cpkj ≥ 0 (see subsection
2.4.2 below) as well as start-up and shutdown costs, cukj , cdkj ≥ 0 (see 2.4.3):
XX
min cpkj + cukj + cdkj .
j∈J k∈K

2.4.2 Production Costs (constraint 2.22 in subsection 2.7)

In [3] a convex production cost function is assumed and approximated by piecewise affine linear
functions. However, since the efficiency of a power unit usually increases with its production level,
we assume a concave production cost function. Fortunately, the increase in efficiency is typically
quite small, allowing a nearly affine linear approximation

∀ j ∈ J, k ∈ K : cpkj = A · L · pkj + B · L · vjk ,

with parameters A and B. A can be interpreted as the variable cost, incurred for every additional
MWh of production, and B as the fixed cost, incurred for running the unit for one hour.
2 THE UNIT COMMITMENT PROBLEM 7

In our model we further split the parameters A and B in two parts, one part which is attributed to
the needed fuel, and is thus dependent on the current fuel price, while the other part is attributed to
operation and maintenance, and thus fixed. We expect the fuel price to be given as the parameter
FCkFj , where Fj 1 denotes the fuel type used by unit j.
Therefore, we replace A and B by their parts, the time and fuel-type dependend FAj and
independend FBj fuel needs and the time and fuel-type dependend PAj and independend PBj
costs. This leads to
 
∀ j ∈ J, k ∈ K : cpkj = FAj · FCkFj + PAj L · pkj
 
+ FBj · FCkFj + PBj L · vjk .

Note that due to the affine, non-linear cost function, we do not have a constant production efficiency
as modeled in many other papers; instead, the modeled production efficiency increases with the
production and is concave, which is typical for thermal units.

2.4.3 Startup and Shutdown Costs (constraints 2.23 and 2.24 in subsection 2.7)

Every unit incurs costs when starting up or shutting down, the former increasing with the offline
time (e.g. for a thermal unit, the start-up costs are partially attributed to the need for reheating it),
the latter constant.
We expect the start-up cost of unit j after an offline time of t periods to be given as CUtj and
the shutdown costs to be given as a single constant CDj . The start-up and shutdown costs can then
be modeled as

∀ j ∈ J, k ∈ [2..T ] : cdkj ≥ CDj (vjk−1 − vjk )

t
!
X
∀ j ∈ J, k ∈ K, t ∈ [k − 1] : cukj ≥ CUtj vjk − vjk−n .
n=1

Pt
Here, the term vjk − n=1 vjk−n is 1 exactly if unit j starts up in period k and was offline in
periods [(k − 1)..(k − t)], otherwise it is less or equal to 0. Thus, the latter constraint is equivalent
to

∀ j ∈ J, k ∈ K, t ∈ {t0 ∈ [k − 1] : ∀ i ∈ [t0 ] : vjk−i = 0} : cukj ≥ CUtj .

Furthermore, since CUtj is increasing with t, this again is equivalent to

∀ j ∈ J, k ∈ K, t = max{t0 ∈ [k − 1] : ∀ i ∈ [t0 ] : vjk−i = 0} : cukj ≥ CUtj .

Since we minimize the costs, cukj = CUtj for start-ups in an optimal solution, as intended.
1F may not contain special characters and spaces, see subsection 4.3
j
2 THE UNIT COMMITMENT PROBLEM 8

2.4.4 Thinning Out the Startup Cost Function

The number of constraints needed in [3] to describe the start-up cost function is quite substantial.
While all the other constraints amount to about 12|J||K|, the model needs about |J||K|2 constraints
to model the start-up cost function as discussed in the last sections. Of course it is not necessary to
model the start-up cost function for |K| cooldown periods, but only for the first periods, when the
cost changes significally. Still, for typical thermal units the relevant timespan is about 2 to 3 days,
leading to between 48|J||K| and 72|J||K| constraints if we assume hourly periods, for example.
If we want to reduce the number of constraints further, we have to use a single step for a whole
group of cooldown periods, which amounts to assigning the same start-up cost to each of these
periods. As an example, let us have a look at a typical start-up cost function with a fixed cost of
about 70% and a variable cost of about 30%. It is possible to reduce the number of steps from 71
to 9, while maintaining a relative error of less than 5%:

100%
Original startup costs (71 steps)
Thinned out startup costs (9 steps)

70%

0% t

Figure 1: Approximation of the start-up cost function with a tolerance of 5%

Thus, two questions are to be answered: How should we group the periods together to obtain a
relative error as small as possible, and which start-up cost should we assign to each of these group?
The second question is answered easily: If we are given a group of periods with cooldown
times ta to tb , it is straight-forward to show that the minimal relative error

CUtjb − CUtja
bestError(CUtja , CUtjb ) =
CUtjb + CUtja
is obtained with the step value
2 · CUtja · CUtjb
bestStep(CUtja , CUtjb ) = ,
CUtjb + CUtja

and 0, if CUtja = CUtjb = 0.

2 THE UNIT COMMITMENT PROBLEM 9

To answer the first question: let STARTUP_TOL be a parameter giving the maximal error tolerance.
Then the periods may be grouped iteratively by starting the first group at t = 1, expanding the
group as long as its relative error is less than the tolerance, and continuing with the next group on
the remaining t’s. In pseudo-code, for given j:

Algorithm 1: startupCostThinning
1 ta , tb ← 1 // Start and end index of current group
2 while ta ≤ |K| do
// Expand the group as long as the next relative error is
less than STARTUP_TOL
3 while tb + 1 ≤ |K| ∧ bestError(CUtja , CUtjb +1 ) < STARTUP_TOL do
4 tb ← tb + 1
5 CUtja ← bestStep(CUtja , CUtjb ) // Calculate optimal step value
6 ∀ t ∈ [(ta +1)..tb ] : CUtj ← 0 // Delete other step values
7 ta , tb ← tb + 1 // Continue with next group

It can be proven that the number of groups produced by this algorithm is the minimal number of
groups needed to fulfill the tolerance requirements.

2.5 Dealing with Infeasibilities

The main cause of infeasibilities is an excessive demand or reserve capacity,
X
∃k: Pj,k < Dj + Rj .
j∈J

In most cases the demand not satisfiable by the units within the model is covered by energy
sources which are not part of the model (e.g. by power imports from external markets or backup
units). Therefore it is advisable to soften the demand and reserve constraints by measuring the
underproduction and underreserve and applying a penalty to it. The softened demand and reserve
constraints thus are:
X
∀k∈K : (pkj − ckj ) + P−k = Dk
j∈J
X
∀k∈K : (pkj − pkj + ckj ) + R−
k
≥ Rk
j∈J

Given the two parameters UPP (penalty factor for underproduction) and URP (penalty factor for
underreserve), we replace the objective function by
XX X
min cpkj + cukj + cdkj + UPP · P−k + URP · R−
k
.
j∈J k∈K k∈K
2 THE UNIT COMMITMENT PROBLEM 10

The second most common cause is a demand profile too volatile to be satisfied by the units, or
in other words, units too slow to follow the demand profile. This infeasibility also originates
from energy sources not present in the model. Although the softening against underproduction
is enough to prevent infeasibilities, one may want to penalize overproduction less (in real world,
overproduction may for example be dumped on external markets). So we obtain
X
∀k∈K : (pkj − ckj ) + P−k − P+k = Dk
j∈J
X
∀k∈K : (pkj − pkj + ckj ) + R−
k
≥ Rk
j∈J

and, given the parameter OPP (penalty factor for overproduction),

XX X
min cpkj + cukj + cdkj + UPP · P−k + URP · R−
k
+ OPP · P+k .
j∈J k∈K k∈K

There are a few other inconsistencies in the input data which could lead to infeasibilities or
unwanted behavior, namely

• simultaneous positive initial down- and uptime:

∃ j ∈ J : IUTj > 0 and IDTj > 0,

• impossible production and ramping limits:

∃ j ∈ J : Pj > Pj or Pj > SUj or Pj > SDj ,

• decreasing start-up costs:

∃ j ∈ J, k ∈ K : CUkj > CUk+1
j , and

• non fulfillable storage constraints:

∃ j ∈ J : SIFj > Pj
∃ j ∈ J : SEj ∈
/ [0, 1]
∃ j ∈ J : SIj > SCj or SFj > SCj
∃ j ∈ J : SIj + L · T · (−SEj Pj + SIFj ) < SFj or SIj + L · T · (−Pj + SIFj ) > SFj

All those are errors in the input data, which have to be avoided for meaningful optimization.
2 THE UNIT COMMITMENT PROBLEM 11

2.6 Parameters summary

2.6.1 General parameters

Name Domain Scale unit Description

L R+ h Period length
UPP R+ cost/MWh Penalty for underproduction
URP R+ cost/MWh Penalty for underreserve
OPP R+ cost/MWh Penalty for overproduction
STARTUP_TOL [0,1] Startup cost error tolerance

2.6.2 Parameters for each period

Name Domain Scale unit Description

Dk R+ MW Power demand
Rk R+ MW Needed reserve capacity
FCkf R+ cost/MWh Cost of a MWhs worth of a fuel type

2.6.3 Parameters for each power unit

Name Domain Scale unit Description

UTj [T ] periods Minimal uptime after starting up
DTj [T ] periods Minimal downtime after shutting down
IUTj [0..T ] periods Initial uptime (unit is online in periods 1, . . . , IUTj )
IDTj [0..T ] periods Initial downtime (unit is offline in periods 1, . . . , IDTj )
Pj R MW Minimal production in online state
Pj R+ MW Maximal production in online state
RUj R+ MW/h Maximal production increase per hour, in online state
RDj R+ MW/h Maximal production decrease per hour, in online state
SUj R+ MW Maximal production increase at start-up
SDj R+ MW Maximal production decrease on shutdown
SCj R+ MWh Storage capacity (Zero for most units)
SEj [0, 1] MWh/MWh Storage efficiency (Zero for most units)
SIFj R+ MW Storage inflow (Zero for most units)
SIj R+ MWh Initial storage fill
SFj R+ MWh Final storage fill
Fj fuel Used fuel type 1
FAj R+ MWh/MWh Fuel need increase per MW of production
FBj R+ MWh/h Fixed fuel need in online state
PAj R+ cost/MWh Cost increase per MWh of production
PBj R+ cost/h Fixed cost in online state
CUtj R+ cost Startup cost after t offline periods
CDj R+ cost Shutdown cost

1F may not contain special characters and spaces, see subsection 4.3
j
2 THE UNIT COMMITMENT PROBLEM 12

2.7 Model summary

Objective function (discussion in 2.4.1 and 2.5, implementation in 3.7):
XX X
min cpkj + cukj + cdkj + UPP · P−k + URP · R−
k
+ OPP · P+k (2.1)
j∈J k∈K k∈K

Physical unit variables (discussion in 2.2.1 and 2.2.5, implementation in 3.2):

∀ j ∈ J, k ∈ K : vjk ∈ {0, 1} (2.2)

∀ j ∈ J, k ∈ K : pkj , pkj ≥ 0 (2.3)
∀ j ∈ J, k ∈ K : skj , ckj ≥0 (2.4)

Cost variables (discussion in 2.4, implementation in 3.2)

∀ j ∈ J, k ∈ K : cpkj , cdkj , cukj ≥ 0 (2.5)

Period variables (discussion in 2.3 and 2.5, implementation in 3.2):

∀k∈K : P−k , P+k , R−

k
≥0 (2.6)

Minimal up- and downtime (discussion in 2.2.2, implementation in 3.3):

∀ j ∈ J, k ∈ [IUTj ] : vjk = 1 (2.7)

∀ j ∈ J, k ∈ [IDTj ] : vjk = 0 (2.8)
∀ j ∈ J, k ∈ [(IUTj + 2)..T ], i ∈ [UTj − 1] ∩ [T − k] : vjk+i ≥ vjk − vjk−1 (2.9)
∀ j ∈ J, k ∈ [(IDTj + 2)..T ], i ∈ [DTj − 1] ∩ [T − k] : vjk+i ≤ 1 − (vjk − vjk−1 ) (2.10)

Minimal and maximal production (discussion in 2.2.3, implementation in 3.3):

∀ j ∈ J, k ∈ K : Pj vjk ≤ pkj ≤ pkj ≤ Pj vjk (2.11)

Ramping constraints (discussion in 2.2.4, implementation in 3.3):

∀ j ∈ J, k ∈ [2..T ] : pkj ≤ pk−1

j + L · RUj vjk−1 + SUj (1 − vjk−1 ) (2.12)
− min{SUj , Pj + L · RUj } · (1 − vjk )
∀ j ∈ J, k ∈ [2..T ] : pkj ≥ pk−1
j − L · RDj vjk − SDj (1 − vjk ) (2.13)
+ min{SDj , Pj + L · RDj } · (1 − vjk−1 )
∀ j ∈ J, k ∈ [T − 1] : pkj ≤ Pj vjk+1 + SDj (vjk − vjk+1 ) (2.14)
2 THE UNIT COMMITMENT PROBLEM 13

Storage constraints (discussion in 2.2.5, implementation in 3.3):

∀ j ∈ J, k ∈ K : skj ≤ SCj (2.15)

∀ j ∈ J, k ∈ K : ckj ≤ max{0, −Pj } (2.16)
∀ j ∈ J, k ∈ [2..T ] : skj = sk−1
j +L· (SEj ck−1
j − pk−1
j + SIFj ) (2.17)
∀j ∈J : s1j = SIj (2.18)
∀j ∈J : SFj = sTj +L· (SEj cTj − pTj + SIFj ) (2.19)

Demand and reserve (discussion in 2.3 and 2.5, implementation in 3.4):

X
∀k∈K : (pkj − ckj ) + P−k − P+k = Dk (2.20)
j∈J
X
∀k∈K : (pkj − pkj + ckj ) + R−
k
≥ Rk (2.21)
j∈J

Production costs (discussion in 2.4.2, implementation in 3.5):

 
∀ j ∈ J, k ∈ K : cpkj = FAj · FCkFj + PAj L · pkj (2.22)
 
+ FBj · FCkFj + PBj L · vjk

Startup and shutdown costs (discussion in 2.4.3 and 2.4.4, implementation in 3.6):

∀ j ∈ J, k ∈ [2..T ] : cdkj ≥ CDj (vjk−1 − vjk ) (2.23)

t
!
X
∀ j ∈ J, k ∈ K, t ∈ [k − 1] : cukj ≥ CUtj vjk − vjk−n (2.24)
n=1
3 IMPLEMENTING THE MODEL 14

3 Implementing the Model

In this section, we show how to implement the model discussed in the last section in Mosel [1],
i.e. how the model composed by parameters, variables and constraints is stated and solved. The
how-to on the surrounding tasks, data input and output, is given in the following sections.
For the whole implementation, we assume the following options to be used:
options noimplicit, explterm, keepassert;

The option noimplicit disallows implicit variable creation, and thus forces the programmer to
declare each variable. Mosel’s automatic line termination requires us to wrap multi-line formulas
after an operator, therefore we deactivate it with explterm. We check the input data with assert(),
which by default is only active in debug mode; keepassert enables it in all modes.

3.1 Parameters (subsection 2.6)

The fitting datatype is clear from the parameters table: real parameters are represented as real,
integer parameters as integer and fuel types as string. The index sets of the parameters are either
J, K, the fuel types or a combination of two of these sets. For performance reasons, the order of the
index sets should be the same as in the forall loops defined over these sets. Here, this means that
the set K always comes last.
Based on these design decisions, we declare the parameters as
declarations
J: set of integer; ! Set of units
K = 1..T; ! Set of periods
Fuels: set of string; ! Available fuel types (index set for CF)
UPP, URP, OPP: real; ! Penalties [cost/MWh]

D, R: array(K) of real; ! Demand and reserve [MW]

FC: dynamic array(Fuels) of
array(K) of real; ! Cost of a fuel type [cost/MWh]

UT, DT: array(J) of integer; ! Minimal up-/downtime [periods]

IUT, IDT: array(J) of integer; ! Initial up-/downtime [periods]
P_min, P_max: array(J) of integer;! Minimal and maximal production [MW]
RU, RD: array(J) of integer; ! Maximal upwards/downwards ramping [MW/h]
SU, SD: array(J) of integer; ! Startup/shutdown ramping [MW]
SC: array(J) of real; ! Storage capacity [MWh]
SE: array(J) of real; ! Storage efficiency [MWh/MWh]
SIF: array(J) of real; ! Storage inflow [MW]
SI, SF: array(J) of real; ! Initial and final storage fill [MWh]
F: array(J) of string; ! Fuel type of a unit [Fuels]
FA, PA: array(J) of real; ! Variable prod. costs [MW/MW],[cost/MW]
FB, PB: array(J) of real; ! Fixed production costs [MW],[cost]
CU: dynamic array(J, range) of real;! Startup costs [cost]
CD: array(J) of real; ! Shutdown costs [cost]
end-declarations
3 IMPLEMENTING THE MODEL 15

The parameters T (number of periods) and L (period length) are declared in the data input section
(4), since they need to be known for selecting the data to be read.
The size of CU’s second index set range is unknown, since we do not know in advance how
many different start-up costs are available for each unit. We therefore declare this array as dynamic.
FC has to be a dynamic array too, since the set of fuels is not known in advance. Usually, we
would declare such an array as
FC: dynamic array(Fuels, K) of real;

However, using the more specific declaration

FC: dynamic array(Fuels) of array(K) of real;

has an important advantage: This way we may pass individual columns of the matrix, the subarrays
FC(f) to the SQL functions, see subsection 4.3. Conveniently, the elements of FC can still be
accessed as FC(f, k).

3.2 Variables
Again, the index sets of the variables are determined from the used indices. In Mosel, the usual
datatype for variables is mpvar (mathematical programming decision variable). The only exception
is the objective function, which as a linear function is declared as linctr (linear constraint, also
used for linear functions).
As the number of storage units is usually quite small, we declare the storage and consumption
variable arrays as dynamic and create the individual variables only for them, thus saving variables.
declarations
p: array(J, K) of mpvar; ! Production [MW]
p_max: array(J, K) of mpvar; ! Maximal possible production [MW]
v: array(J, K) of mpvar; ! On-Off state [0/1]
s: dynamic array(J, K) of mpvar; ! Storage [MWh]
c: dynamic array(J, K) of mpvar; ! Consumption [MW]
cp: array(J, K) of mpvar; ! Production costs [cost]
cu, cd: array(J, K) of mpvar; ! Startup/shutdown costs [cost]
overallCosts: linctr; ! Overall costs [cost]
p_under, p_over: array(K) of mpvar; ! Under-/overproduction [MW]
r_under: array(K) of mpvar; ! Underreserve [MW]
end-declarations

! Create storage and consumption variables for storage units

forall(j in J, k in K | P_min(j) < 0) do
create(s(j, k));
create(c(j, k));
end-do

By default, these variables are constrained to take positive real values. This fits all variables except
the on-off-state v, which should be binary:
forall(j in J, k in K) v(j, k) is_binary;
3 IMPLEMENTING THE MODEL 16

3.3 Physical Constraints of the Units (constraints 2.7 to 2.14)

The Mosel equivalent of the minimal up- and downtime constraints (2.7 to 2.10) is canonical:
forall(j in J, k in 1..IUT(j)) v(j, k) = 1;
forall(j in J, k in 1..IDT(j)) v(j, k) = 0;

for the initial up- and downtime, and

forall(j in J, k in IUT(j)+2 .. T, i in 1..minlist(UT(j)-1, T-k))
v(j, k + i) >= v(j, k) - v(j, k-1);
forall(j in J, k in IDT(j)+2 .. T, i in 1..minlist(DT(j)-1, T-k))
v(j, k + i) <= 1 - (v(j, k-1) - v(j, k));

for the interperiod minimal up- and downtime.

The jointed constraints for minimal and maximal production (2.11) have to be separated:
forall(j in J, k in K) do
P_min(j) * v(j, k) <= p(j, k);
p(j, k) <= p_max(j, k);
p_max(j, k) <= P_max(j) * v(j, k);
end-do

The ramping limits (2.12 to 2.14) are straight-forward to implement:

forall(j in J, k in 2..T) do
p_max(j, k) <= p(j, k-1) + L * RU(j) * v(j, k-1) + SU(j) * (1 - v(j, k-1))
- minlist(SU(j), P_min(j) + L * RU(j)) * (1 - v(j, k));
p(j, k) >= p(j, k-1) - L * RD(j) * v(j, k) - SD(j) * (1 - v(j, k))
+ minlist(SD(j), P_min(j) + L * RD(j)) * (1 - v(j, k-1));
end-do
forall(j in J, k in 1..T-1) do
p_max(j, k) <= P_max(j) * v(j, k+1) + SD(j) * (v(j, k) - v(j, k+1));
end-do

Same as with the storage variables, the storage constraints (2.15 to 2.19) are implemented only
for storage units.
forall(j in J | P_min(j) < 0) do
forall(k in K) do
s(j, k) <= SC(j);
c(j, k) <= maxlist(0, -P_min(j));
end-do
forall(k in 2..T) do
s(j, k) = s(j, k-1) + L * (SE(j)*c(j, k-1) - p(j, k-1) + SIF(j));
end-do
SI(j) = s(j, 1);
SF(j) = s(j, T) + L * (SE(j)*c(j, T) - p(j, T) + SIF(j));
end-do
3 IMPLEMENTING THE MODEL 17

3.4 Power Grid Constraints (constraints 2.20 and 2.21)

Since the consumption variables are only created for storage units, we use exists to sum only the
existing variables.
forall(k in K) do
sum(j in J) p(j, k) - sum(j in J | exists(c(j, k))) c(j, k)
+ p_under(k) - p_over(k) = D(k);
sum(j in J) (p_max(j, k) - p(j, k)) - sum(j in J | exists(c(j, k))) c(j, k)
+ r_under(k) >= R(k);
end-do

3.5 Production Cost Constraints (constraint 2.22)

The model could be reduced by replacing every use of cpkj by the right-hand side of constraints 2.22,
and thus removing the variables cpkj completely. Fortunately, the Xpress Optimizer automatically
applies this reduction at the presolve stage and enables us to use the canonical implementation for
better readability:
forall(j in J, k in K) cp(j, k) = (FA(j) * FC(F(j), k) + PA(j)) * p(j, k)
+ (FB(j) * FC(F(j), k) + PB(j)) * v(j, k);

3.6 Startup and Shutdown Costs (constraints 2.23 and 2.24)

The start-up and shutdown constraints can be implemented in Mosel as
forall(j in J, k in 2..T) cd(j, k) >= CD(j) * (v(j, k-1) - v(j, k));
forall(j in J, k in K, t in 1..k-1 | exists(CU(j, t)))
cu(j, k) >= CU(j, t) * (v(j, k) - sum(n in 1..t) v(j, k - n));

The exists operator is again used to enumerate only relevant indices. Especially the number
of constraints depending on t may be further reduced from considering only a subset (compare
subsection 2.4.4), accepting a loss in accuracy of the start-up costs.

3.6.1 Thinning Out the Startup Cost Function

The pseudo-code of the start-up cost thinning can be translated one-to-one into Mosel code:
function bestError(CU_j_ta: real, CU_j_tb: real): real
if(CU_j_ta = 0 and CU_j_tb = 0) then
returned := 0;
else
returned := abs(CU_j_ta - CU_j_tb)/(CU_j_ta + CU_j_tb);
end-if
end-function
3 IMPLEMENTING THE MODEL 18

function bestStep(CU_j_ta: real, CU_j_tb: real): real

if(CU_j_ta = 0 and CU_j_tb = 0) then
returned := 0;
else
returned := 2 * CU_j_ta * CU_j_tb / (CU_j_ta + CU_j_tb);
end-if
end-function

procedure startupCostThinning(j: integer)

declarations
t_a, t_b: integer;
end-declarations
! Start the first group
t_a := 1;
t_b := 1;
repeat
! Expand the group as long as the next
! relative error is less than STARTUP_TOL
while(exists(CU(j,t_b + 1)) and
bestError(CU(j,t_a),CU(j,t_b+1)) < STARTUP_TOL) do
t_b := t_b+1;
end-do
! Calculate optimal step value
CU(j, t_a) := bestStep(CU(j, t_a), CU(j, t_b));
! Delete other step values
forall(t in t_a+1..t_b) do
delcell(CU(j, t));
end-do
! Continue with next group
t_a := t_b + 1;
t_b := t_b + 1;
until(not exists(CU(j, t_a)));
end-procedure

Of course, the thinning has to be applied to every unit j ∈ J:

forall(j in J) do
startupCostThinning(j);
end-do

3.7 Objective Function (constraint 2.1)

First, we have to define the cost function
overallCosts := sum(j in J, k in K) (cp(j, k) + cd(j, k) + cu(j, k))
+ sum(k in K) (UPP*p_under(k) + URP*r_under(k) + OPP*p_over(k));

and then call minimize:

minimize(overallCosts);
3 IMPLEMENTING THE MODEL 19

3.8 Dealing with Infeasibilities (subsection 2.5)

The penalties have already been incorporated in the power grid constraints and the objective
function. To detect errors in the input data, we use Mosel’s assert procedure, which stops the
program if a condition is not met.
assert(and(j in J) IUT(j) in 0..T, "Initial uptime out of range!");
assert(and(j in J) IDT(j) in 0..T, "Initial uptime out of range!");
assert(and(j in J) IUT(j)*IDT(j) = 0, "Simultaneous initial down- and uptime!");
assert(and(j in J) UT(j) in 1..T, "Initial uptime out of range!");
assert(and(j in J) DT(j) in 1..T, "Initial downtime out of range!");
assert(and(j in J) P_min(j) <= P_max(j), "Impossible production limits!");
assert(and(j in J) P_min(j) <= SU(j), "Some unit is not able to start up!");
assert(and(j in J) P_min(j) <= SD(j), "Some unit is not able to shutdown!");
declarations
lastCU: real;
end-declarations
forall(j in J) do
lastCU := 0;
forall(k in K | exists(CU(j, k))) do
assert(lastCU <= CU(j, k),
"The start-up costs are not monotonically increasing!");
lastCU := CU(j, k);
end-do
end-do
assert(and(j in J) SIF(j) <= P_max(j), "Storage inflow leads to overcapacity!");
assert(and(j in J) (0 <= SE(j) and SE(j) <= 1), "Invalid storage efficiency!");
assert(and(j in J) SI(j) <= SC(j), "Invalid initial storage fill!");
assert(and(j in J) SF(j) <= SC(j), "Invalid final storage fill!");
assert(and(j in J) SI(j) + L*T*(SE(j)*maxlist(0, -P_min(j)) + SIF(j)) >= SF(j),
"Some storage constraints are not fulfillable!");
assert(and(j in J ) SI(j) + L*T*(-P_max(j) + SIF(j)) <= SF(j),
"Some storage constraints are not fulfillable!");

As noted at the beginning of this section (see 3), we need to use the option keepassert to activate
assert outside of debug mode.
4 DATA INPUT 20

4 Data Input
4.1 General Parameters
We use the two parameters UNITS_SOURCE and PERIODS_SOURCE to store the sources for the unit
parameters and the period data. Depending on whether one wants to access a database or a
spreadsheet, these parameters contain the name of a database table or a spreadsheet.
Once we know where to get the period data, we have to know a start time START and the number
of periods T. Since the datatype datetime can not be used for parameters, START will be a string,
formatted with the standard SQL format YYYY-MM-DD HH:MM:SS (example: June 30th, 2007
at 7:12 PM is written as 2007-06-30 19:12:00). The parameter PERIOD_LENGTH should be given in
seconds.
Summarising, the used parameters are:
parameters
UNITS_SOURCE = "Units.xls"; ! Source of the unit parameters
PERIODS_SOURCE = "Periods.xls"; ! Source of the period data
DSN = "Server=?;Database=?;UID=?;PWD=?;"; ! Example DSN
START = "2009-05-11 00:00:00"; ! Start of the modeled timespan
T = 168; ! Number of modeled periods
L = 1; ! Period length [h]
STARTUP_TOL = 0.05; ! Tolerance for modeling of start-up cost
end-parameters

Important: parameters have to be declared with a default value to define their data type.

4.2 Unit Parameters

The parameters describing the power units are too many to be read from the command line, so we
have to read them from a different data source. Mosel supports database access, including Excel
spreadsheets, over ODBC.
Now, databases are superior to spreadsheets in performance, reliability and scalability, making
them a good choice for a Mosel model used in day-to-day operations. Excel spreadsheets on the
other hand are easy to create and modify, making them a good choice for the development and
experimental stage.
Fortunately, the Mosel initializations block hides most of the differences between real
databases and Excel spreadsheets, allowing us to switch between them with little effort. For details
on how to use the initializations block, please refer to the FICO Whitepaper “Using ODBC
and other database interfaces with Mosel”, which should reside in the Xpress installation directory
at .../XpressMP/docs/mosel/mosel_odbc/moselodbc.pdf [4].
We expect two tables in our database,

• one table with the unit parameters with index set J (all except start-up costs),

• one table with the start-up costs with index sets J and K.
4 DATA INPUT 21

The name of the first table should be given as the parameter UNITS_SOURCE, whereas the second
table should have the same name with an appended “_CU” (in reference to the start-up costs name
CUtj ). The names of the index and of the columns should be the same as in our model.
The Excel spreadsheet is set up in the same way, except that we store each set of units in their
own spreadsheet file. Thus the UNITS_SOURCE parameter now denotes the spreadsheet file, while
the two sheets are always called “Units” and “Units_CU”.
Finally, we need the ODBC data source name (DSN), also called “Connection String”, which
differs with the brand of database server used. We expect the DSN to be given as DSN. Connection
strings can be found in the database’s manual. A collection of connection strings for popular
databases is also available at www.connectionstrings.com.
Now, we can setup the parameters needed by initializations:
declarations
unitsDSN: string; ! ODBC data source name
unitsTable: string; ! Table containing the unit parameter
unitsCUTable: string; ! Table containing the start-up costs
end-declarations
if(isExcelDocument(UNITS_SOURCE)) then
unitsDSN := "mmodbc.excel:skiph;" + UNITS_SOURCE;
unitsTable := "Units";
else
unitsDSN := "mmodbc.odbc:" + DSN;
unitsTable := UNITS_SOURCE;
end-if
unitsCUTable := unitsTable + "_CU";

The function isExcelDocument used here is listed in Appendix A, and works by checking the file
extension.
In the initializations block, we read all unit parameters with index set J from the first table
initializations from unitsDSN
[UT,DT,IUT,IDT,P_min,P_max,RU,RD,SU,SD,SC,SE,SIF,SI,SF,F,FA,FB,PA,PB,CD]
as unitsTable +
"(j,UT,DT,IUT,IDT,P_min,P_max,RU,RD,SU,SD,SC,SE,SIF,SI,SF,F,FA,FB,PA,PB,CD)";

and the start-up costs from the second table:

CU as unitsCUTable + "(j,k,CU)";
end-initializations

Mosel automatically uses the columns j and k as the unit and period indices. Once the units have
been read, we are able to define the Fuels set:
Fuels := union(j in J) {F(j)};
finalize(Fuels);
4 DATA INPUT 22

4.3 Period Data

Since we only want to get data for the periods in our timespan, not the full dataset, and we (possibly)
do not know all used fuel types, the number of arrays to be read must be variable, which is not
possible within one initializations. An elegant solution to this problem is the direct access to
the database using Mosel’s SQL functions.
The period data is stored in a single table which holds the demand D, the reserve R and the
fuel costs. For every used fuel, we expect the apposite fuel cost column to be called FC_f, where
f stands for the actual name of the fuel type. Therefore, the fuel name may not contain special
characters and spaces.
To access the database, we have to setup a connection first, so the ODBC data source name is
needed again:
if(isExcelDocument(PERIODS_SOURCE)) then
SQLconnect("DSN=Excel Files;HDR=Yes;DBQ="+expandpath(PERIODS_SOURCE));
assert(getparam("SQLsuccess"), ’SQLconnect failed!’);
else
SQLconnect(DSN);
assert(getparam("SQLsuccess"), ’SQLconnect failed!’);
end-if

Note that expandpath expands the path of PERIODS_SOURCE, since the ODBC driver’s working
directory might be different from ours.
To send an SQL query to the database now, we have to build it up first. The following variables
are used for this:
declarations
periodsIndex: text; ! Index column
reindexedPeriods: text; ! Period source with converted index
periodsColumns: list of text; ! Columns to read
periodsArrays: list of array(K) of real;! Arrays to fill
periodsQuery: text; ! Assembled SQL query
end-declarations

The first column to be read is the period index k. How to derive this index depends on the indexing
of the periods in the database. For our model, we expect the database periods to be indexed by the
column t with data type TIMESTAMP. We also assume the database periods to have the same length
as the model periods, L; data with different period lengths needs to be interpolated beforehand.
Given all these informations, we can derive our index k as
periodsIndex := "{fn FLOOR(((t-{ts ’"+START+"’})*24+0.1)/"+L+")} + 1";

• {ts ’"+START+"’} is the start time START as a timestamp

• (t-{ts ’"+START+"’}) is the elapsed time since START in days
• ((t-{ts ’"+START+"’})*24+0.1) is the elapsed time in hours
L
• (((t-{ts ’"+START+"’})*24+0.1)/"+L+") is the elapsed time in periods + 10
• {fn FLOOR(((t-{ts ’"+START+"’})*24+0.1)/"+L+")} is the period index ∈ [0..T − 1]
• {fn FLOOR(((t-{ts ’"+START+"’})*24+0.1)/"+L+")} + 1 is our index k ∈ [T ]
4 DATA INPUT 23

In a fully-fledged database, it would be possible to assign the name k to this converted index, and
to use it by this name. In Excel, this assignment has to be done in a subquery. At the same time,
we have to specify the source of the periods:
reindexedPeriods := "SELECT "+periodsIndex+" as k, *";
if(isExcelDocument(PERIODS_SOURCE)) then
reindexedPeriods += " FROM Periods";
else
reindexedPeriods += " FROM " + PERIODS_SOURCE;
end-if

The subquery reindexedPeriods now represents the periods table, with a converted and renamed
index column k.
With the index settled, we can select our first data arrays for reading: the demand D and the
reserve R.
periodsColumns += ["D", "R"];
periodsArrays += [D, R];

Here, the last statement actually does not copy the D and R arrays. Since they are complex datatypes,
Mosel just copies their references to the list. So, when initialising the arrays of this list, we actually
initialize our original arrays.
The procedure for the fuel costs is similar, but since the FC array is dynamic, we first have to
create the cost subarray for each used fuel type:
forall(f in Fuels) do
create(FC(f));
periodsColumns += ["FC_" + f];
periodsArrays += [FC(f)];
end-do

Finally, we can assemble the SQL query and specify the needed periods:
periodsQuery := "SELECT k, " + join(periodsColumns, ", ");
periodsQuery += " FROM ("+reindexedPeriods+")";
periodsQuery += " WHERE 1 <= k AND k <= " + T;

The join function joins the elements of periodsColumns, its implementation can be found in
Appendix A.
In conclusion, an assembled query for example may equate to
SELECT k, D, R, FC_gas, FC_oil, FC_coal
FROM (
SELECT {fn FLOOR(((t-{ts ’2009-05-11 00:00:00’})*24+0.1)/1)} + 1 as k, *
FROM Periods
)
WHERE 1 <= k AND k <= 84

We now initialize our parameters using SQLexecute:

SQLexecute(periodsQuery, periodsArrays);
assert(getparam("SQLsuccess"), ’SQLexecute failed!’);
4 DATA INPUT 24

After the initialization, we have to close the connection to the database server:
SQLdisconnect;

When working with big Excel spreadsheets, this approach is a good alternative to multiple
initializations blocks. Excel reopens the spreadsheet for every initializations block, which
can be quite time consuming. In contrast, when using SQL functions, Excel opens the spreadsheet
just once at SQLconnect.
5 DATA OUTPUT 25

5 Data Output
Same as with the data input sources, the destination file is given as a parameter:
parameters
RESULTS_DEST = "Results.xls"; ! Destination of the results
end-parameters

5.1 Electricity price

Once the cost function has been minimized within the model, we can derive a basic estimate of the
electricity price by observing:
• The price is set the highest accepted bid of any unit.

• In a market with perfect competition, a unit will bid in at marginal costs.

• The marginal cost of a unit is FAj · FCkFj + PAj .

This can be calculated in Mosel as:
declarations
price: array(K) of real;
end-declarations
forall(k in K) do
price(k) := max(j in J | v(j, k).sol >= 0.5) (FA(j) * FC(F(j), k) + PA(j));
end-do

5.2 Postprocessing
The maximal production and ramping constraints just impose an upper limit on the maximal
possible production. Thus, p_max may actually be lower than the real maximal possible production.
We can fix this by deriving the maximal possible production from the production variables and the
operational state:
declarations
exact_p_max: array(J, K) of real;
end-declarations
forall(j in J, k in K) do
exact_p_max(j, k) := P_max(j).sol * v(j, k).sol;
if(k > 1) then
exact_p_max(j, k) := minlist(exact_p_max(j, k),
p(j, k-1).sol + RU(j) * v(j, k-1).sol
+ SU(j) * (1 - v(j, k-1).sol) + P_max(j) * (1 - v(j, k).sol));
end-if
if(k < T) then
exact_p_max(j, k) := minlist(exact_p_max(j, k),
P_max(j) * v(j, k+1).sol + SD(j) * (v(j, k).sol - v(j, k+1).sol));
end-if
end-do
5 DATA OUTPUT 26

5.3 Data Output to Excel

We want to output the optimizer variables for

1. the operational state vjk ,

2. the production variables pkj , pkj ,

3. the price pricej and

4. the costs cpkj , CUkj , CDj,k .

We will output each variable to its own sheet (inside the same file). This allows us to

• change the sizes of J and K and

• add new variables to the model

without having to change the position of the existing variables. Since there are usually more periods
than units, and Excel supports more rows than columns2 , we associate the periods with rows and
the units with columns. Together with the actual data, we have to output the index sets J and K.
Since the period indices are relative to START and L, we will also output the timestamp of each
period.
Now, it would be convenient to implement the output in a separate function, which is to be
called for each variable. Unfortunately, one initializations block per function call would be
needed, reopening the Excel spreadsheet each time. Depending on the size of the spreadsheet,
this may be quite time-consuming. Therefore we use a single initializations block, but still
perform the preparation of each variable in a separate function.
We need two functions: one for variables with index set K, and one for variables index sets J
and K. They both should return a two-dimensional array of text, which subsequently will be
written to an Excel sheet:
declarations
PeriodsSheet = dynamic array(rows: range, singleCol: range) of text;
UnitsPeriodsSheet = dynamic array(rows, cols: range) of text;
end-declarations

function createSheet(data: array(K) of real, name: string): PeriodsSheet

...
end-function

function createSheet(data: array(J, K) of real): UnitsPeriodsSheet

...
end-function

2 Excel 2007 or higher: 220 rows vs. 214 columns; prior to Excel 2007: 216 rows vs. 28 columns or less
5 DATA OUTPUT 27

In the body of the createSheet function for variables with index set K, we need to:

1. Fill the first two columns with the index and the starting time of each period:
forall(k in K) do
returned(k + 1, 1) := text(k);
returned(k + 1, 2) := text(datetime(START) + (k-1) * L * 3600);
end-do

2. Fill the first row with the given variable name:

returned(1, 3) := name;

3. Fill the space between period indices and titles with the actual data:
forall(k in K) do
returned(k + 1, 3) := text(data(k));
end-do

In the function for variables with index sets J and K, the time series for all units J should be
written to adjacent columns. Since the index set J is not necessarily of the form {1, . . . , |J|}, we
first construct a mapping from J to {1, . . . , |J|}:
declarations
mapping: dynamic array(J) of integer;
end-declarations
forall(j in J) do
mapping(j) := getsize(mapping) + 1;
end-do

The remainder of the function is similiar to the previous function for index set K. We need to:

1. Fill the first two columns with the index and the starting time of each period:
forall(k in K) do
returned(k + 1, 1) := text(k);
returned(k + 1, 2) := text(datetime(START) + (k-1) * L * 3600);
end-do

2. Fill the first row with the indices of the units,

forall(j in J) do
returned(1, mapping(j)) := text(j);
end-do

3. Fill the space between period indices and titles with the actual data:
forall(j in J, k in K) do
returned(k + 1, mapping(j)) := text(data(j, k));
end-do
5 DATA OUTPUT 28

These functions only take real-valued time series; For convenience, we implement them for
mpvar-valued time series as well:
function createSheet(data: array(K) of mpvar, name: string): PeriodsSheet
returned := createSheet(array(k in K) data(k).sol, name);
end-function

function createSheet(data: array(J, K) of mpvar): UnitsPeriodsSheet

returned := createSheet(array(j in J, k in K) data(j, k).sol);
end-function

To perform the output to Excel, we still need to specify the region to which each variable needs
to be written. For this, we will use the cellRange function, implemented in our ExcelFunctions
package (Appendix A), which generates Excel ranges like [Sheet$A1:Q10].
initializations to "mmodbc.excel:noindex;" + RESULTS_DEST
evaluation of createSheet(v)
as cellRange("v", 1, 1, getsize(J)+2, getsize(K)+1);
evaluation of createSheet(p)
as cellRange("p", 1, 1, getsize(J)+2, getsize(K)+1);
evaluation of createSheet(exact_p_max)
as cellRange("p_max", 1, 1, getsize(J)+2, getsize(K)+1);
evaluation of createSheet(price, "Price")
as cellRange("price", 1, 1, 3, getsize(K)+1);
evaluation of createSheet(s)
as cellRange("s", 1, 1, getsize(J)+2, getsize(K)+1);
evaluation of createSheet(c)
as cellRange("c", 1, 1, getsize(J)+2, getsize(K)+1);
evaluation of createSheet(cp)
as cellRange("cp", 1, 1, getsize(J)+2, getsize(K)+1);
evaluation of createSheet(cu)
as cellRange("cu", 1, 1, getsize(J)+2, getsize(K)+1);
evaluation of createSheet(cd)
as cellRange("cd", 1, 1, getsize(J)+2, getsize(K)+1);
end-initializations
6 REMARKS 29

6 Remarks
This whitepaper presents a detailled model for the Unit Commitment problem. We have given its
practical implementation and discussed how to reduce the number of start-up cost constraints and
how to analyze the cause of infeasibilities. Finally, we have shown how to manage data input and
output, and have given a basic estimate of the electricity price.
One major drawback of the implemented model is the lacking consideration of market power.
To remedy this, usually a so-called “uplift” function is added to the model’s electricity price, which
is derived by statistically analyzing the disagreement between the outcome of this model and the
real price on historical data.
Of course, it would be better to consider the market power directly as part of the model.
Different approaches in achieving this are discussed in [5]. The two most popular ones are game
theoretical approaches:

• the Cournot equilibrium and

• the Supply Function equilibrium.

They differ in the modelling of the strategies of the power generating companies. While in a
Cournot model a company may only decide on it’s power output, a company’s strategy in a Supply
Function equilibrium is described by a function mapping the market price to its power supply. The
greater freedom in the choice of a strategy, however, comes at the expense of a higher computational
effort, which in turn forces the use of less detailled models.
For a comparison of the practical performance of Cournot and Supply Function models on the
German power market, we refer to [6].
The second drawback is the neglection of the underlying power grid. The production schedule
resulting from our model may not be feasible on a real power grid, i.e. the power grid may not be
able to transport the electricity from the producing units to the consumers. While this was typically
not a fundamental issue so far, it becomes more and more important due to the increasing energy
production from renewable, more volatile resources causing bottlenecks in the power exchange
beetween different regions and thus necessitate an explicit modelling of the power grid.
Also, most of today’s power markets are connected to one or more neighboring markets through
interconnectors. The interconnectors are used to transport energy from a market with higher price
to a market with lower price, thus diminishing the price difference. The effective price change is
typically small due to the small capacity of the interconnectors, but may not be negligible.
REFERENCES 30

References
[1] FICO Xpress-Mosel. http://www.fico.com/en/Products/DMTools/xpress-
overview/Pages/Xpress-Mosel.aspx.

[2] FICO Xpress Optimization Suite. http://www.fico.com/en/products/

fico-xpress-optimization-suite/.

[3] M. Carrión and J.M. Arroyo. A computationally efficient mixed-integer linear formulation for
the thermal unit commitment problem. Power Systems, IEEE Transactions on, 21(3):1371–
1378, 2006.

[4] S. Heipcke. Using odbc and other database interfaces with mosel. 2009. http://www.
fico.com/en/wp-content/secure_upload/Xpress_moselodbc.pdf.

[5] M. Ventosa, A. Baillo, A. Ramos, and M. Rivier. Electricity market modeling trends. Energy
policy, 33(7):897–913, 2005.

[6] B. Willems, I. Rumiantseva, and H. Weigt. Cournot versus supply functions: What does the
data tell us? Energy Economics, 31(1):38–47, 2009.
A EXCEL FUNCTIONS 31

A Excel Functions
The Excel functions are implemented as part of the ExcelFunctions package. Therefore, the
functions to be used in other packages or the main model have to marked as public.

A.1 Function isExcelDocument

The function isExcelDocument checks if the given data source is an Excel file. This is done by
checking the file extension:
public function isExcelDocument(source: string): boolean
declarations
S = getsize(source);
end-declarations
returned := strtolower(copytext(text(source), S-4+1, S)) = ’.xls’ or
strtolower(copytext(text(source), S-5+1, S)) = ’.xlsx’;
end-function

Unfortunately, we also have to implement the strtolower function for ourselfs, which converts a
string to lower case:
public function strtolower(str: text): text
returned := str;
forall(i in 1..getsize(returned)) do
if(65(! = A !) <= getchar(returned, i) and
getchar(returned, i) <= 90(! = Z !)) then
setchar(returned, i, getchar(returned, i) + (97-65));
end-if
end-do
end-function

A.2 Referencing cells in Excel

Excel supports two notations for referencing cells,

• the A1 notation (ex: C4, R7, Z80) and

• the R1C1 notation (ex: R4C3, R7C18, R80C26).

While the R1C1 notation is easier to use programmatically, the A1 notation may be more familiar
to end users. Therefore, we offer the boolean useA1notation which switches from R1C1 to A1
notation:
public declarations
useA1notation: boolean; ! If true, the A1 notation is used
end-declarations

Depending on this setting, the coordinates of a cell can be determined with a call to cellCoords:
A EXCEL FUNCTIONS 32

public function cellCoords(column: integer, row: integer): string

if(useA1notation) then
returned := columnToA1(column) + row;
else
returned := "R" + row + "C" + column;
end-if
end-function

(→ columnToA1 is described in the next subsection)

Cell ranges are described by the coordinates of their top-left and bottom-right cells, divided by a
colon (ex. C4:R7, R4C3:R7C18). They can be determined with a call to cellRangeCoords:
public function cellRangeCoords(column: integer, row: integer,
width: integer, height: integer): string
returned := cellCoords(column, row) + ":"
+ cellCoords(column + width - 1, row + height - 1);
end-function

When referencing Excel cells through ODBC, it is usual to include the cell’s sheet name too. This
is done by prepending the sheet name to the coordinates, divided by a dollar sign, and by enclosing
the reference in square brackets. The two following functions generate such references to a cell or
respectively to a cell range:
public function cell(sheet: string, column: integer, row: integer): string
returned := "[" + sheet + "$" + cellCoords(column, row) + "]";
end-function
public function cellRange(sheet: string, column: integer, row: integer,
width: integer, height: integer): string
returned := "["+sheet+"$"+cellRangeCoords(column, row, width, height) + "]";
end-function

A.2.1 Column index conversions in the A1 notation

The column numbering scheme in the A1 notation is a bit unusual and needs special conversion
functions.
Conversion from integer to A1 column index (1 → A, 2 → B, ...)
public function columnToA1(index: integer): string
declarations
c: text; ! Column name as text (needs mmsystem)
end-declarations
while(index > 0) do
c := " " + c;
setchar(c, 1, 65 (! = A !) + ( (index-1) mod 26));
index := (index-1) div 26;
end-do
returned := string(c);
end-function

Conversion from A1 column index to integer (A → 1, B → 2, ...)

A EXCEL FUNCTIONS 33

public function a1ToColumn(column: string): integer

forall(i in 1..getsize(column)) do
returned := returned * 26 + getchar(column, i) - 65 (! = A !) + 1;
end-do
end-function

A.3 Function Join

Joins the specified separator string between each element of the specified set of strings, yielding a
single joined string. Example:
join(["A", "D", "B", "C"], ", ") -> "A, D, B, C"
join(["A"], ", ") -> "A"
join({"A", "D", "B", "C"}, ", ") -> "A, B, C, D"

Implementations for lists and sets of string:

public function join(strings: list of text, separator: string): text
forall(s in strings) do
returned := returned + s + separator;
end-do
returned -= separator;
end-function

public function join(strings: set of text, separator: string): text

forall(s in strings) do
returned := returned + s + separator;
end-do
returned -= separator;
end-function