# Energy-Saving Generation Dispatch Using Minimum Cost Flow

Authors: Zhan’an Zhang; Xingguo Cai
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

This study uses a minimum cost flow method to solve a dispatch problem in order to minimize the consumption of coal in the dispatching of a thermal power system. Low-carbon generation dispatching is also considered here since the scheduling results are consistent with energy-saving generation dispatch. Additionally, this study employs minimum coal consumption as an objective function in considering the output constraints, load balance constraints, line loss, ramp rate limits, spinning reverse needs, prohibited operating zone requirements, security constraints, and other common constraints. The minimum cost flow problem, considering the loss of network flow, is known as a generalized network flow problem, which can be expressed as a quadratic programming problem in mathematics. Accordingly, the objective function was solved by LINGO11, which was used to calculate a network in a single time; a continuous period dispatch plan was obtained by accumulating each period network flow together. This analysis proves it feasible to solve a minimal cost flow problem with LINGO11. Theoretical analysis and numerical results prove the correctness and effectiveness of the proposed method.


---

## Body

## 1. Introduction

There seems to be rather compelling evidence that global warming is an issue that we seriously need to be concerned about today [1, 2]. Carbon dioxide accounts for 80% of the greenhouse effect, and rising carbon dioxide levels are the main cause of global warming [3, 4]. China pledged to reduce its carbon intensity by 40–45% by 2020 based on 2005 levels. By the end of 2012, 71.5% of China’s generation capacity was from thermal power, of which 92.5% was the product of coal-fired generation [5]. Coal-fired electricity consumes about 50% of China’s coal production, and the CO2 emissions from power generation account for 40% of the total CO2 emissions in China. Accordingly, China’s power industry has implemented energy-saving generation dispatching (ESGD) and low-carbon generation dispatching (LCGD).ESGD is one of the most important problems in power system operations requiring load demand at minimum total fuel cost while accounting for various unit and system constraints. The ESGD model is an optimization problem that considers linear and nonlinear characteristics, including power balance constraints, generation limit constraints, node voltage constraints, ramp rate limits, spinning reverse needs, prohibited operating zone requirements, and security constraints, among others [6]. In this paper, low-carbon generation dispatching is also considered, since the scheduling results are consistent with energy-saving generation dispatching. In practical scheduling applications, a daily scheduling period is generally divided into 24 or more intervals; therefore, dispatching of each period can be solved as a static optimization problem.Many mathematical techniques have been developed and applied to dispatch problem such as linear programming [7], interior-point method [8], Lagrangian relaxation algorithm [9], quadratic programming [10] and other traditional algorithms. These algorithms essentially need some problem simplification such that the problem is linear or convex. Thus, a true global minimum cannot be guaranteed [11]. The dynamic programming method [12] has also been successfully used in solving the dispatch problems; however, this method may result in “curse of dimensionality.” More recently, the metaheuristic algorithms such as particle swarm optimization (PSO) [13], genetic algorithm (GA) [14], simulated annealing (SA) [15], evolutionary programming (EP) [16], and ant colony optimization (ACO) [17], have also been considered in the context of economic dispatch. However, the biggest problem that metaheuristic algorithms faced is that the computational efficiency is rather low. The related parameters are not easy to set up and the computational time is long. It is not easy to implement when applied to large electrical power systems. Their highly heuristic nature usually leads to suboptimal solutions. Differential evolution (DE) is a stochastic search based method [18], which can present a simple structure, fast convergence speed, and robustness. However, DE fast convergence might lead the direction of the search toward a local optimal and premature solution. Essentially, the economic dispatch problem is a large scale nonlinear programming problem. In pursuit of the optimal solution for economic dispatch, various hybrid methods have been investigated and implemented [5, 19–22]. Unfortunately, these hybrid algorithms normally take lengthy calculation time when compared with the mathematical optimization methods. Moreover, previous algorithms rarely considered network structures, we use a network flow method here to make up for the deficiency, because the network flow method can well retain the topology of the network.Networks provide a useful way for modeling power system problems and are used extensively in power system dispatching [23, 24]. As an important network problem, ESGD problems can be formulated and solved as minimum cost flow problems when the cost is a quadratic function of the power, which is a nonlinear, minimum cost flow problem.The main objective of this study is to introduce convex quadratic programming to solve the ESGD problem, since coal consumption and network losses are all convex functions of the power flow through a network. In order to do that this study employed LINGO11 to solve quadratic programming problems accounting for linear and nonlinear equality constraints and inequality constraints. Accordingly, a minimum cost flow algorithm was used to solve the ESGD model and calculate a network within a single moment. From that process, a continuous period of ESGD planning was obtained by accumulating period network flow results. This process confirmed that the minimal cost flow method was successful for solving the ESGD problem and, therefore, has value for these types of applications.


## 2. Problem Formulation

### 2.1. The Mathematical Model of ESGD Problem

The ESGD problem determines the optimal schedule of the available generating units to simultaneously minimize the generation cost rate and meet the load demand of a power system while meeting various equality and inequality constraints. This mathematical model can be described as follows:(1)
m
i
n
∑
t
=
1
T
∑
i
=
1
N
F
i
t
P
i
t
=
∑
t
=
1
T
∑
i
=
1
N
α
i
P
i
t
2
+
β
i
P
i
t
+
γ
i
,where F
i
t is the cost function of the ith generator; α
i, β
i, and γ
i are the cost coefficients of the ith generator; P
i
t is the power of the ith generator at t time; and N represents the number of generators committed to the operating system.


### 2.2. The Mathematical Model of Low-Carbon Generation Dispatching

Achieving the lowest carbon emission is the target of low-carbon generation dispatching. This mathematical model can be described as follows:(2)
m
i
n
∑
t
=
1
T
∑
i
=
1
N
d
i
t
P
i
t
,where d
i
t
(
P
i
t
) is the electrical-carbon characteristic function. This formula represents the CO2 emissions when the output of unit i is P
i
t at t time (t/h), which can be expressed as (3)
d
t
=
e
q
η
P
t
=
2.77038
F
P
t
,where e is the CO2 emission coefficient of the fuel used in a power source; the standard coal emission factor is 2.77, which means that 2.77 kg CO2 can be discharged for every 1 kg of standard coal burnt; q is the calorific value of unit fuel, which is 8.14 kWh/kg of standard coal; and η is power generation efficiency, which can be expressed as(4)
η
=
3600
P
t
29308
F
P
t
×
1
0
3
×
100
%
,where 3600 is the electric heating value (kJ/kWh) and 29308 is the calorific value of standard coal (kJ/kg).If the coal consumption function and electrical-carbon characteristic satisfies the following relationship(5)
F
1
P
1
t
d
1
P
1
t
=
⋯
=
F
i
P
1
t
d
i
P
1
t
=
λ
,
d
i
P
1
t
≠
0
,
F
i
P
1
t
=
0
,
d
i
P
1
t
=
0
,where λ is a constant, then we can get (6)
m
i
n
∑
t
=
1
T
∑
i
=
1
N
F
i
t
P
i
t
=
λ
m
i
n
∑
t
=
1
T
∑
i
=
1
N
d
i
t
P
i
t
;then the 2 kinds of scheduling results are consistent [25].


### 2.3. Constraints

The objective function needs to satisfy the following constraints.(1) Power Constraints. The power constraints include generator output, transformer capacity, and line transmission limits. Consider(7)
P
i
m
i
n
≤
P
i
,
t
≤
P
i
m
a
x
i
=
1,2
,
…
,
N
,where P
i
m
i
n is the lower and P
i
m
a
x is the upper output limit of unit i, respectively.(2) Node Voltage Constraints. Consider(8)
U
i
m
i
n
≤
U
i
,
t
≤
U
i
m
a
x
i
=
1,2
,
…
,
N
,where U
i
,
t is the voltage of node i at t time and U
i
m
i
n is the lower and U
i
m
a
x is the upper limit of node i, respectively.(3) Power Balance Constraints. Consider(9)
∑
i
=
1
N
P
i
,
t
=
P
D
t
+
P
L
t
i
=
1,2
,
…
,
N
,where P
D is the total load demand and P
L is the transmission network losses, which is a function of unit power outputs that can be represented using the B coefficients:(10)
P
L
=
∑
i
=
1
N
∑
j
=
1
N
P
i
B
i
j
P
j
+
∑
i
=
1
N
P
i
B
0
i
+
B
00
,where B
i
j is the loss coefficient square matrix; B
0
i is the loss coefficient vector; and B
00 is the loss coefficient constant [26].(4) Operation Ramp Rate Limits. The power output of a practical generator cannot be adjusted instantaneously without limits. The operating range for all online units is restricted by their ramp rate limits during each dispatch period. Therefore, the dispatch output of a generator should be limited by the constraints of up and down ramp rates [27], which are given as follows:(11)
P
i
,
t
-
1
-
P
i
,
t
≤
D
R
i
i
f
g
e
n
e
r
a
t
i
o
n
d
e
c
r
e
a
s
e
s
P
i
,
t
-
P
i
,
t
-
1
≤
U
R
i
i
f
g
e
n
e
r
a
t
i
o
n
i
n
c
r
e
a
s
e
s
,where D
R
i and U
R
i are the ramp-down and ramp-up rate limits of the ith thermal unit, respectively [28]. If the unit ramp rate limits are considered, the real power operating limits are modified as follows:(12)
m
a
x
P
i
m
i
n
,
P
i
,
t
-
1
-
D
R
i
≤
P
i
,
t
≤
m
i
n
P
i
m
a
x
,
P
i
,
t
-
1
+
U
R
i
.(5) Spinning Reserve. The added spinning reserve factor must be considered to prevent a sudden large load to the system or a failure in a certain large unit requirement. This condition can be explained as follows:(13)
∑
i
=
1
N
P
i
m
a
x
×
U
i
t
≥
P
D
t
+
R
t
,where R
t is the spinning reserve in the tth hour and U
i
t is the ON and OFF status of the ith conventional unit at the t period (U
i
t
=
0 represents OFF status and U
i
t
=
1 represents ON status) [6].(6) Prohibited Operating Zone. The prohibited operating zone is the range of prohibited output power resulting from the physical limitations of machine components, steam valves, vibration in the shaft bearing, and other conditions that can cause discontinuity in the electrical energy cost curve. Therefore, some units must be considered as prohibited zones in practical operation. The feasible operating zones of thermal units can be described as follows [29]: (14)
P
i
m
i
n
≤
P
i
≤
P
i
,
1
l
,
P
i
,
j
-
1
u
≤
P
i
≤
P
i
,
j
l
,
P
i
,
z
i
u
≤
P
i
≤
P
i
m
a
x
,
j
=
2,3
,
…
,
z
i
,where z
i is the number of prohibited zones in the ith generator curve; j is the index of prohibited zone of the ith generator; and P
i
,
j
u and P
i
,
j
l are the upper and lower limits of the jth prohibited zone of unit i, respectively.(7) Security Constraints. For secure operation, the transmission line loading (S
i) is restricted by its upper limit as follows:(15)
S
i
≤
S
i
m
a
x
i
=
1,2
,
…
,
N
L
,where N
L is the total number of lines [30].


## 2.1. The Mathematical Model of ESGD Problem

The ESGD problem determines the optimal schedule of the available generating units to simultaneously minimize the generation cost rate and meet the load demand of a power system while meeting various equality and inequality constraints. This mathematical model can be described as follows:(1)
m
i
n
∑
t
=
1
T
∑
i
=
1
N
F
i
t
P
i
t
=
∑
t
=
1
T
∑
i
=
1
N
α
i
P
i
t
2
+
β
i
P
i
t
+
γ
i
,where F
i
t is the cost function of the ith generator; α
i, β
i, and γ
i are the cost coefficients of the ith generator; P
i
t is the power of the ith generator at t time; and N represents the number of generators committed to the operating system.


## 2.2. The Mathematical Model of Low-Carbon Generation Dispatching

Achieving the lowest carbon emission is the target of low-carbon generation dispatching. This mathematical model can be described as follows:(2)
m
i
n
∑
t
=
1
T
∑
i
=
1
N
d
i
t
P
i
t
,where d
i
t
(
P
i
t
) is the electrical-carbon characteristic function. This formula represents the CO2 emissions when the output of unit i is P
i
t at t time (t/h), which can be expressed as (3)
d
t
=
e
q
η
P
t
=
2.77038
F
P
t
,where e is the CO2 emission coefficient of the fuel used in a power source; the standard coal emission factor is 2.77, which means that 2.77 kg CO2 can be discharged for every 1 kg of standard coal burnt; q is the calorific value of unit fuel, which is 8.14 kWh/kg of standard coal; and η is power generation efficiency, which can be expressed as(4)
η
=
3600
P
t
29308
F
P
t
×
1
0
3
×
100
%
,where 3600 is the electric heating value (kJ/kWh) and 29308 is the calorific value of standard coal (kJ/kg).If the coal consumption function and electrical-carbon characteristic satisfies the following relationship(5)
F
1
P
1
t
d
1
P
1
t
=
⋯
=
F
i
P
1
t
d
i
P
1
t
=
λ
,
d
i
P
1
t
≠
0
,
F
i
P
1
t
=
0
,
d
i
P
1
t
=
0
,where λ is a constant, then we can get (6)
m
i
n
∑
t
=
1
T
∑
i
=
1
N
F
i
t
P
i
t
=
λ
m
i
n
∑
t
=
1
T
∑
i
=
1
N
d
i
t
P
i
t
;then the 2 kinds of scheduling results are consistent [25].


## 2.3. Constraints

The objective function needs to satisfy the following constraints.(1) Power Constraints. The power constraints include generator output, transformer capacity, and line transmission limits. Consider(7)
P
i
m
i
n
≤
P
i
,
t
≤
P
i
m
a
x
i
=
1,2
,
…
,
N
,where P
i
m
i
n is the lower and P
i
m
a
x is the upper output limit of unit i, respectively.(2) Node Voltage Constraints. Consider(8)
U
i
m
i
n
≤
U
i
,
t
≤
U
i
m
a
x
i
=
1,2
,
…
,
N
,where U
i
,
t is the voltage of node i at t time and U
i
m
i
n is the lower and U
i
m
a
x is the upper limit of node i, respectively.(3) Power Balance Constraints. Consider(9)
∑
i
=
1
N
P
i
,
t
=
P
D
t
+
P
L
t
i
=
1,2
,
…
,
N
,where P
D is the total load demand and P
L is the transmission network losses, which is a function of unit power outputs that can be represented using the B coefficients:(10)
P
L
=
∑
i
=
1
N
∑
j
=
1
N
P
i
B
i
j
P
j
+
∑
i
=
1
N
P
i
B
0
i
+
B
00
,where B
i
j is the loss coefficient square matrix; B
0
i is the loss coefficient vector; and B
00 is the loss coefficient constant [26].(4) Operation Ramp Rate Limits. The power output of a practical generator cannot be adjusted instantaneously without limits. The operating range for all online units is restricted by their ramp rate limits during each dispatch period. Therefore, the dispatch output of a generator should be limited by the constraints of up and down ramp rates [27], which are given as follows:(11)
P
i
,
t
-
1
-
P
i
,
t
≤
D
R
i
i
f
g
e
n
e
r
a
t
i
o
n
d
e
c
r
e
a
s
e
s
P
i
,
t
-
P
i
,
t
-
1
≤
U
R
i
i
f
g
e
n
e
r
a
t
i
o
n
i
n
c
r
e
a
s
e
s
,where D
R
i and U
R
i are the ramp-down and ramp-up rate limits of the ith thermal unit, respectively [28]. If the unit ramp rate limits are considered, the real power operating limits are modified as follows:(12)
m
a
x
P
i
m
i
n
,
P
i
,
t
-
1
-
D
R
i
≤
P
i
,
t
≤
m
i
n
P
i
m
a
x
,
P
i
,
t
-
1
+
U
R
i
.(5) Spinning Reserve. The added spinning reserve factor must be considered to prevent a sudden large load to the system or a failure in a certain large unit requirement. This condition can be explained as follows:(13)
∑
i
=
1
N
P
i
m
a
x
×
U
i
t
≥
P
D
t
+
R
t
,where R
t is the spinning reserve in the tth hour and U
i
t is the ON and OFF status of the ith conventional unit at the t period (U
i
t
=
0 represents OFF status and U
i
t
=
1 represents ON status) [6].(6) Prohibited Operating Zone. The prohibited operating zone is the range of prohibited output power resulting from the physical limitations of machine components, steam valves, vibration in the shaft bearing, and other conditions that can cause discontinuity in the electrical energy cost curve. Therefore, some units must be considered as prohibited zones in practical operation. The feasible operating zones of thermal units can be described as follows [29]: (14)
P
i
m
i
n
≤
P
i
≤
P
i
,
1
l
,
P
i
,
j
-
1
u
≤
P
i
≤
P
i
,
j
l
,
P
i
,
z
i
u
≤
P
i
≤
P
i
m
a
x
,
j
=
2,3
,
…
,
z
i
,where z
i is the number of prohibited zones in the ith generator curve; j is the index of prohibited zone of the ith generator; and P
i
,
j
u and P
i
,
j
l are the upper and lower limits of the jth prohibited zone of unit i, respectively.(7) Security Constraints. For secure operation, the transmission line loading (S
i) is restricted by its upper limit as follows:(15)
S
i
≤
S
i
m
a
x
i
=
1,2
,
…
,
N
L
,where N
L is the total number of lines [30].


## 3. The Minimal Cost Flow Method

### 3.1. Network Flow Theory Introduction

Network flow problems can be described with graph theory, where a graph represents a network of nodes and connecting arcs. If a path exists between any two pairs of vertices in a graph, then that graph is a connected graph where each arc has a specific direction.As depicted in Figure1, the common abstraction that models a flow network is a directed graph G
=
(
V
,
E
), where V is the set of vertices and E is the set of edges over these vertices.Figure 1
A network flow diagram.A
r
c
(
i
,
j
) is the edge connecting nodes i and j, which has a predetermined direction, flow f
i
j, the cost of the unit flow c
i
j, and flow rate limit (as shown in Figure 1). Each edge (
i
,
j
) has a flow f
i
j that defines the number of units of the commodity that flows from i to j. An edge also has a capacity limit that constrains the number of units that can flow over that edge. Here, c
i
j denotes the cost per unit flow transported directly from vertex i to vertex j [31]. When no units are flowing over an edge, then the arc is outlined in a dotted line.A special source vertexs
∈
V produces units of a commodity that flow through the edges of the graph to be consumed by a sink vertex r
∈
V known as the receiving node. The rest nodes, for which there are both outflow and inflow arcs, are termed intermediate points.The following criteria must be satisfied for any feasible flow through a network.(1) Capacity Constraint. The flow f
i
j through an edge cannot exceed the capacity of the edge u
i
j,  l
i
j
≤
f
i
j
≤
u
i
j. If an edge (
i
,
j
) does not exist in the network, then u
i
j
=
0.(2) Flow Conservation. Aside from the source vertex s and sink vertex r, each vertex i
∈
V must satisfy the property for which the sum of f
i
j for all edges (
i
,
j
) in E (the flow into i) must equal the sum of f
i
j for all edges (
i
,
j
)
∈
E (the flow out of i). This property ensures that flow is neither produced nor consumed in the network, except at s and r.(3) Skew Symmetry. For consistency, the quantity f
i
j represents the net flow from vertex i to j. This means that it must be the case that f
i
j
=
-
f
j
i, which holds, even if both edges (
i
,
j
) and (
j
,
i
) exist in a directed graph (Figure 1) [32].


### 3.2. The Minimal Cost Flow Method

The minimum cost flow problem determines the minimum total cost under the specified flow while taking into account the arc cost. The objective function can be expressed as follows:(16)
m
i
n
∑
(
i
,
j
)
∈
E
c
i
j
f
i
j
.Ifc
i
j is a constant, the minimum cost flow problem can be solved efficiently since it can be formulated as a linear programming problem.Both the coal consumption and the network losses are convex functions of power flow within a power network. In minimal cost flow problems, the cost of the arc is a convex function of the flow on the arc and that is a convex cost flow problem. Whenc
i
j is a linear function of flow f
i
j, then the objective function changes into a quadratic programming problem. Not only the objective function but also the constraints include quadratic terms of f
i
j. Accordingly, the cost function of unit flow can be expressed as(17)
d
f
P
i
d
P
=
2
α
i
P
i
+
β
i
,where f
(
P
i
) stands for F
(
P
i
) or P
L. Hence, the objective function model can be expressed as follows:(18)
m
i
n
∑
i
=
1
n
α
i
f
i
+
β
i
f
=
m
i
n
∑
i
=
1
n
α
i
f
i
2
+
β
i
f
i
s
.
t
.
∑
i
,
j
∈
E
f
i
j
-
μ
j
i
∑
j
,
i
∈
E
f
j
i
=
b
i
,
∀
i
∈
G
.In this model, the value ofb
(
i
) depends on the nature of node i, where b
(
i
)
>
0 if node i is a supply node, b
(
i
)
<
0 if node i is a demand node, and b
(
i
)
=
0 when node i is a transshipment node [33].In a general minimum cost flow problem, the arc is the conservation of the flow, and the flow entering an arc equals the flow leaving the arc. This assumption is reasonable in many practical application scenarios; however, the power flow will diminish when flowing through the grid due to resistance, which, in this study, is considered an issue associated with generalized flow problems. In generalized flow problems, arcs might consume or generate flow. Iff
i
j units of flow enter an a
r
c
(
i
,
j
), then μ
i
j
f
i
j units arrive at node j, where μ
i
j is a positive multiplier associated with the arc. If 0
≤
μ
i
j
≤
1, the arc is lossy, and, if 1
≤
μ
i
j
≤
∞, the arc is gainy. The problem becomes a general minimal cost flow problem when μ
j
i
=
1 and b
(
i
)
=
0. In this study, all the arcs are loss arcs due to the resistance, which means that (
0
≤
μ
i
j
≤
1
) [34].There is no single standard algorithm that can always be used to solve convex programming problems. There have been many algorithms for solving convex quadratic programming problems, such as the Lemke method, the interior-point method [35, 36], the effective set method [37, 38], and an ellipsoid algorithm, for example, each having its own advantages and disadvantages.In this study, LINGO11, very mature software widely used in various mathematical optimization problems, is used to solve the quadratic programming algorithm problems.


### 3.3. Multiperiod Network Flow Model

A continuous time network flow is constituted by combining each network flow chart for a single period of time. A multiperiod model is a three-dimensional model with each layer corresponding to a single-period model (as shown in Figure2). Each s
i is united to form new s, and, in the same way, to form a new r. A
r
c
(
s
,
s
i
) represents the output and a
r
c
(
r
i
,
r
) represents the load. Considering each period together, a continuous-time power scheduling plan can be obtained. It is merely a schematic diagram, and the practical applications in power system may be more complex than that in Figure 2. For instance, the direction of the arcs in each layer may vary with different operating conditions. In order to simplify the analysis, we assume that the lines and transformers meet their capacity constraints.Figure 2
Multiperiod network flow model.The above structure model (Figure2) represents a typical generalized network flow programming problem. It can be expressed as(19)
m
i
n
C
T
X
,satisfied with(20)
A
X
=
b
,where C represents the m-dimensional arc cost vector; X represents the m-dimensional arc flow; b represents an n-dimensional vector injected nodes; and A represents n
×
m dimensional node-arc incidence matrix [39].


## 3.1. Network Flow Theory Introduction

Network flow problems can be described with graph theory, where a graph represents a network of nodes and connecting arcs. If a path exists between any two pairs of vertices in a graph, then that graph is a connected graph where each arc has a specific direction.As depicted in Figure1, the common abstraction that models a flow network is a directed graph G
=
(
V
,
E
), where V is the set of vertices and E is the set of edges over these vertices.Figure 1
A network flow diagram.A
r
c
(
i
,
j
) is the edge connecting nodes i and j, which has a predetermined direction, flow f
i
j, the cost of the unit flow c
i
j, and flow rate limit (as shown in Figure 1). Each edge (
i
,
j
) has a flow f
i
j that defines the number of units of the commodity that flows from i to j. An edge also has a capacity limit that constrains the number of units that can flow over that edge. Here, c
i
j denotes the cost per unit flow transported directly from vertex i to vertex j [31]. When no units are flowing over an edge, then the arc is outlined in a dotted line.A special source vertexs
∈
V produces units of a commodity that flow through the edges of the graph to be consumed by a sink vertex r
∈
V known as the receiving node. The rest nodes, for which there are both outflow and inflow arcs, are termed intermediate points.The following criteria must be satisfied for any feasible flow through a network.(1) Capacity Constraint. The flow f
i
j through an edge cannot exceed the capacity of the edge u
i
j,  l
i
j
≤
f
i
j
≤
u
i
j. If an edge (
i
,
j
) does not exist in the network, then u
i
j
=
0.(2) Flow Conservation. Aside from the source vertex s and sink vertex r, each vertex i
∈
V must satisfy the property for which the sum of f
i
j for all edges (
i
,
j
) in E (the flow into i) must equal the sum of f
i
j for all edges (
i
,
j
)
∈
E (the flow out of i). This property ensures that flow is neither produced nor consumed in the network, except at s and r.(3) Skew Symmetry. For consistency, the quantity f
i
j represents the net flow from vertex i to j. This means that it must be the case that f
i
j
=
-
f
j
i, which holds, even if both edges (
i
,
j
) and (
j
,
i
) exist in a directed graph (Figure 1) [32].


## 3.2. The Minimal Cost Flow Method

The minimum cost flow problem determines the minimum total cost under the specified flow while taking into account the arc cost. The objective function can be expressed as follows:(16)
m
i
n
∑
(
i
,
j
)
∈
E
c
i
j
f
i
j
.Ifc
i
j is a constant, the minimum cost flow problem can be solved efficiently since it can be formulated as a linear programming problem.Both the coal consumption and the network losses are convex functions of power flow within a power network. In minimal cost flow problems, the cost of the arc is a convex function of the flow on the arc and that is a convex cost flow problem. Whenc
i
j is a linear function of flow f
i
j, then the objective function changes into a quadratic programming problem. Not only the objective function but also the constraints include quadratic terms of f
i
j. Accordingly, the cost function of unit flow can be expressed as(17)
d
f
P
i
d
P
=
2
α
i
P
i
+
β
i
,where f
(
P
i
) stands for F
(
P
i
) or P
L. Hence, the objective function model can be expressed as follows:(18)
m
i
n
∑
i
=
1
n
α
i
f
i
+
β
i
f
=
m
i
n
∑
i
=
1
n
α
i
f
i
2
+
β
i
f
i
s
.
t
.
∑
i
,
j
∈
E
f
i
j
-
μ
j
i
∑
j
,
i
∈
E
f
j
i
=
b
i
,
∀
i
∈
G
.In this model, the value ofb
(
i
) depends on the nature of node i, where b
(
i
)
>
0 if node i is a supply node, b
(
i
)
<
0 if node i is a demand node, and b
(
i
)
=
0 when node i is a transshipment node [33].In a general minimum cost flow problem, the arc is the conservation of the flow, and the flow entering an arc equals the flow leaving the arc. This assumption is reasonable in many practical application scenarios; however, the power flow will diminish when flowing through the grid due to resistance, which, in this study, is considered an issue associated with generalized flow problems. In generalized flow problems, arcs might consume or generate flow. Iff
i
j units of flow enter an a
r
c
(
i
,
j
), then μ
i
j
f
i
j units arrive at node j, where μ
i
j is a positive multiplier associated with the arc. If 0
≤
μ
i
j
≤
1, the arc is lossy, and, if 1
≤
μ
i
j
≤
∞, the arc is gainy. The problem becomes a general minimal cost flow problem when μ
j
i
=
1 and b
(
i
)
=
0. In this study, all the arcs are loss arcs due to the resistance, which means that (
0
≤
μ
i
j
≤
1
) [34].There is no single standard algorithm that can always be used to solve convex programming problems. There have been many algorithms for solving convex quadratic programming problems, such as the Lemke method, the interior-point method [35, 36], the effective set method [37, 38], and an ellipsoid algorithm, for example, each having its own advantages and disadvantages.In this study, LINGO11, very mature software widely used in various mathematical optimization problems, is used to solve the quadratic programming algorithm problems.


## 3.3. Multiperiod Network Flow Model

A continuous time network flow is constituted by combining each network flow chart for a single period of time. A multiperiod model is a three-dimensional model with each layer corresponding to a single-period model (as shown in Figure2). Each s
i is united to form new s, and, in the same way, to form a new r. A
r
c
(
s
,
s
i
) represents the output and a
r
c
(
r
i
,
r
) represents the load. Considering each period together, a continuous-time power scheduling plan can be obtained. It is merely a schematic diagram, and the practical applications in power system may be more complex than that in Figure 2. For instance, the direction of the arcs in each layer may vary with different operating conditions. In order to simplify the analysis, we assume that the lines and transformers meet their capacity constraints.Figure 2
Multiperiod network flow model.The above structure model (Figure2) represents a typical generalized network flow programming problem. It can be expressed as(19)
m
i
n
C
T
X
,satisfied with(20)
A
X
=
b
,where C represents the m-dimensional arc cost vector; X represents the m-dimensional arc flow; b represents an n-dimensional vector injected nodes; and A represents n
×
m dimensional node-arc incidence matrix [39].


## 4. The Solution Process Diagram

The steps in the procedure of ESGD (the low-carbon dispatching calculation process is the same) can be described as follows (as shown in Figure3).Figure 3
Flowchart of solution procedure.Step 1.
Determine the loadP
t at t time according to the load forecasting curve, where T represents the total of time periods.Step 2.
PutP
t as the flow at sink vertex.Step 3.
Solve the solution according to (1), satisfying the constraints (see (7) to (15)).Step 4.
Ift
=
T, then output ESGD scheduling results; otherwise, if t
<
T, then let t
=
t
+
1 and go to Step 2.


## 5. Example Analysis

For this study, we calculated a modified IEEE 6-bus power system, which contains 2 generators, 2 transformers and 5 transmission lines (as shown in Figure4).Figure 4
IEEE 6-bus power system diagram.The parameters in the modified system are shown in Tables1, 2, and 3. In the tables, r
0 is the resistivity of a conductor per unit length, L is the length of lines, R and R
T are the resistance of lines and transformers, respectively, Δ
P
k represents short circuit loss of a transformer, and α
i, β
i, and γ
i refer to the cost coefficients of the ith generator (in Table 3).Table 1
The parameters of transmission lines.
Line
V (kV)
Model
r
0 (Ω/km)
L (km)
R (Ω)
L
1
110
LGJ-185
0.17
71.2
12.1
L
2
110
LGJ-185
0.17
71.2
12.1
L
3
110
LGJ-185
0.17
71.2
12.1
L
4
10
LGJ-120
0.242
10
2.42
L
5
10
LGJ-120
0.242
10
2.42Table 2
The parameters of transformers.
T
V (kV)
Model
Δ
P
k (kw)
R
T (Ω)
T
1
110/11
SFPL1-50000/110
250
1.21
T
2
110/11
SFPL1-50000/110
250
1.21Table 3
The parameters of generators.
G
Capacity (MW)
Output limit
α
i
β
i
γ
i
G
1
100
50 ≤P
G
1 ≤ 100
0.0014
0.25
2.5
G
2
100
50 ≤P
G
2 ≤ 100
0.0018
0.18
5.0The corresponding cost flow network diagram is shown in Figure5, in which the directed diagram of the system includes 8 vertexes and 12 arcs.Figure 5
Network diagram.In order to consider network losses, the net loss unit is turned into a standard coal consumption unit; that is, 1 MW = 0.123 t/h. In this model, the equivalent calorific value of 1 kg standard coal is 7000 kcal. The total-cost optimization model can be expressed as follows:(21)
m
i
n
F
x
=
(
0.25
x
01
+
0.0028
x
01
2
)
h
h
h
h
h
i
+
(
0.18
x
02
+
0.0036
x
02
2
)
h
h
h
h
h
i
+
0.0002
×
0.123
×
(
x
43
2
+
x
65
2
)
h
h
h
h
h
i
+
0.0004
×
0.123
h
h
h
h
h
i
×
(
x
25
2
+
x
23
2
)
+
0.002
×
0.123
h
h
h
h
h
i
×
(
x
14
2
+
x
46
2
+
x
16
2
)
s
.
t
.
x
01
-
x
14
-
x
16
=
0
x
02
-
x
25
-
x
23
=
0
x
23
-
0.0002
x
23
2
+
x
43
-
0.0001
x
43
2
-
x
37
=
0
x
14
-
0.001
x
14
2
-
x
46
-
x
43
=
0
x
65
-
0.0001
x
65
2
+
x
25
-
0.0002
x
25
2
-
x
57
=
0
x
46
-
0.001
x
46
2
+
x
16
-
0.001
x
16
2
-
x
65
-
x
67
=
0
x
37
+
x
57
+
x
67
=
P
D
3
+
P
D
5
+
P
D
6
50
≤
x
01
≤
100
,
50
≤
x
02
≤
100
,
0
≤
x
14
≤
100
,
0
≤
x
46
≤
100
,
0
≤
x
16
≤
100
,
0
≤
x
43
≤
50
,
0
≤
x
25
≤
100
,
0
≤
x
65
≤
50
,
0
≤
x
23
≤
100
,
0
≤
x
67
≤
100
,
0
≤
x
57
≤
100
,
0
≤
x
37
≤
100
,where x refers to the flow in the network; then the objective function is obviously a convex cost flow problem and a quadratic programming method can be used to solve it. When the load changes, as shown in Figure 8, the differing distribution of the minimum cost flow can be concluded and, therefore, a continuous dispatch result can be obtained.In order to simplify the analysis, we assumed that all units run within a safe operation area. The ramp rate limit of 100 MW unit is 2 MW/min, so the ramp rate limits can meet the requirements created when the time interval is an hour.We also assumed that the reactive power can be compensated locally and the system does not transmit reactive power. The voltage of nodes was kept within their specified ratings, and the reserve capacity of the grid was assumed to meet operational requirements. Additionally, the test was performed on an Intel (R) Pentium (R) CPU 2.13 GHz 2.0 GB RAM with LINGO11, and the average iteration times were 58. The maximum flow that can pass the network and the corresponding minimal cost distribution is shown in Figure6. The minimum flow distribution of the minimum cost flow is shown in Figure 7.Figure 6
The minimal cost flow under maximum flow.Figure 7
The minimal cost flow distribution under minimum flow.Figure 8
Load and coal consumption curve.The load and corresponding minimal coal consumption are listed in Table4 and shown in Figure 8, respectively. Where T is the time period (h), P is the power of load (MW) and F is the coal consumption (t/h).Table 4
Load and coal consumption.
T (h)
1
2
3
4
5
6
7
8
9
10
11
12
P (MW)
110
100
110
120
130
140
150
160
170
160
150
140
F (t/h)
44.58
38.74
44.58
50.80
57.39
64.36
71.71
79.45
87.58
79.45
71.71
64.36
T (h)
13
14
15
16
17
18
19
20
21
22
23
24
P (MW)
130
140
150
160
170
180
190
180
170
150
130
120
F (t/h)
57.39
64.36
71.71
79.45
87.58
96.10
105.0
96.10
87.58
71.71
57.39
50.80Electricity and minimal coal consumption accumulation are listed in Table5 and shown in Figure 10, respectively. S refers to the accumulation of coal consumption (t) and W refers to accumulation of electricity (MWh).Table 5
Electricity and coal consumption accumulation.
T (h)
1
2
3
4
5
6
7
8
9
10
11
12
S (t)
44.6
83.3
127.9
178.7
236.1
300.5
372.2
451.6
539.2
618.6
690.4
754.7
W (MWh)
110
210
320
440
570
710
860
1020
1190
1350
1500
1640
T (h)
13
14
15
16
17
18
19
20
21
22
23
24
S (t)
812.1
876.5
948.2
1027
1115
1211
1316
1412
1500
1571
1629
1680
W (MWh)
1770
1910
2060
2220
2390
2570
2760
2940
3110
3260
3390
3510It can be seen that the accumulated electricity will gradually increase with the increase of time. The same trend occurs with the accumulated minimal coal consumption. Unit scheduling results are shown in Figure9.Figure 9
Outputs of generators.Figure 10
Electricity and coal consumption accumulation curve.According to the electrical-carbon conversion relations, 2.77 tons of CO2 can be emitted for 1 ton of standard coal consumed. In addition, the minimum CO2 emissions and the accumulation curve can also be obtained (as shown in Figures 11 and 12).Figure 11
CO2 emissions curve.Figure 12
CO2 emissions accumulation curve.The value of CO2 emissions and its accumulation are listed in Table 6, where C
e and C
a refer to CO2 emissions and its accumulation value, respectively.Table 6
The value of CO2 emissions and CO2 emissions accumulation.
T(h)
1
2
3
4
5
6
7
8
9
10
11
12
C
e (t)
123.5
107.3
123.5
140.7
159.0
178.3
198.6
220.1
242.6
220.1
198.6
178.3
C
a (t)
123.5
230.8
354.3
495.0
654.0
832.3
1031
1251
1494
1714
1912
2091
T (h)
13
14
15
16
17
18
19
20
21
22
23
24
C
e (t)
159.0
178.3
198.6
220.1
242.6
266.2
290.9
266.2
242.6
198.6
159.0
140.7
C
a (t)
2249
2428
2626
2846
3089
3355
3646
3912
4155
4353
4513
4653We can see that the emissions of CO2 change with the power load by the same regularity, which reflects that the objectives of low-carbon generation dispatching and energy-saving generation dispatching are consistent.In order to support our findings, we also solved this same problem by invoking the CPLEX solver in the MATLAB toolbox YALMIP [40]. The extended solver of LINGO11 includes Branch-and-Bound solver, Global solver, and Multistart solver, while the solver in CPLEX is mainly based on interior-point method. If these two different methods come to a same conclusion, then the result concluded is credible. The result by CPLEX produced the same solution garnered using LINGO11, which proves the correctness of the method used in this study.


## 6. Conclusions

The minimal cost flow algorithm can be used to dispatch the power system and can reflect the system network topology, making it fairly easy to consider system constraints. The method is simple, rapid, and clear in terms basic concepts. The example analysis proves that the method is feasible and practical.In order to simplify the analysis, this study only considers the energy-saving scheduling of thermal units, without considering the effects of hydropower, wind power, and other renewable energy sources: although they do not consume coal and discharge CO2, they do have an impact on load-flow distribution and network losses in power systems. Accordingly, the next step will be to initiate research to study the influences of renewable energy sources on ESGD in power systems.



---

*Source: 562462-2015-10-26.xml*