Piecewise Linear Approximation of Generators Cost

Functions Using Max-Affine Functions

Hamed Ahmadi Ali Moshref
José R. Martí BBA
School of Electrical and Computer Engineering Vancouver, BC, Canada
University of British Columbia Email: ali.moshref@bba.ca
Vancouver, BC, Canada
Email: hameda@ece.ubc.ca, jrms@ece.ubc.ca

Abstract—Nonlinear functions are often encountered in power operation is of crucial importance.

system optimizations. In this paper, an effective piecewise linear Most of the optimization problems in power systems are
(PWL) approximation technique is introduced which shows basically nonlinear, thus nonlinear programming (NLP) tech-
promising performance in linearizing the nonlinear functions.
This method uses a series of linear functions, called max-affine niques have been widely applied to these problems. As an
functions, to linearize a multivariate function over a bounded example for continuous optimization, the optimal power flow
domain. The important advantage of this method is its ability to problem considering AC constraints is an NLP problem,
decide on the size of the subspaces, which other methods are not for which a variety of methods have been proposed in the
capable of. It is also shown that using the PWL approximation, literature, e.g. Interior Point Method [1] and Trust-Region
significant efficiency is achievable in computation burden of most
power system optimizations, such as unit commitment. Method [2]. However, due to the nonconvexity of the original
problem, most of the methods would probably get trapped in a
N OMENCLATURE local minimum. As an example of combinatorial optimization,
B̄ Network susceptance matrix. the unit commitment problem is nonconvex due to the binary
D̄ Matrix product of the susceptance and node- variables associated with the generators’ status (on/off). Hav-
incident matrices. ing a quadratic cost function and nonlinear constraints, one has
a, b, c Quadratic cost function coefficients. to deal with a mixed-integer nonlinear programming (MINLP)
C SDn Shutdown cost. problem, which is NP-hard, and up to now, there is no efficient
C SUp Startup cost. method for solving large-scale MINLP problems.
CT System total cost. Various versions of the unit commitment problem have
Nb , Ng Sets of buses and generators, respectively. been formulated in the literature, focusing on linearizing the
Nl , Nt Sets of lines and time horizon, respectively. constraints as well as the objective function, e.g. [3], [4].
P Generator active power. In [3], the quadratic cost functions are linearized within the
P max , P min Upper/lower limits on generator active power. generator’s output power limits which leads to 2s + 2 new
Pd Active power demand. constraints and s + 1 new variables (s is the number of
PiL Active power flow limit of Line i. sections assumed for linearization). Beside increasing the size
P SDn Generator minimum power limit at shutdown. of the problem, it is not clear how to choose the intervals
P SUp Generator maximum power limit at startup. and, therefore, it might not lead to the best possible piecewise
RDn Generator power ramp-down limit. linear (PWL) approximation of the function. In [4], the convex
RUp Generator power ramp-up limit. envelope of the quadratic function is obtained and using the
s Number of PWL partitions. perspective cut, the linear parts are added to the problem. The
SDn, SUp Auxiliary variables for shutdown/startup cost. problem with this method is that a dynamic constraint has
u Generator status (0:’Off’, 1:’On’). to be added to the original problem which slows down the
δ Bus voltage angle. solution process. In addition, a linearization technique is used
in [5] to iteratively solve the MILP problem by replacing the
I. I NTRODUCTION objective by its linear approximation. The method is known as
The increasing demand on real-time operation of power the Kelley’s theorem on the cutting plane method for convex
systems has led to numerous research studies on enhancing the programs. However, many iterations and cuts are required to
power system analysis methods. Due to large-scale problems reach to the solution and the number of required cuts is not
that need to be solved for real-size power systems, special predictable beforehand.
numerical methods have been developed and utilized by power In a more recent work, a mathematical approach is em-
system experts. Amongst most computationally expensive ployed to find the PWL approximation of the quadratic cost
problems in power system analysis, optimization of system functions [6]. Using this method, the length of the intervals can
be selected optimally so that a good PWL approximation is in the approximation. In order to clarify this and without loss
obtained. Although this is believed to be a tighter approxima- of generality, let the dimension of x be one. In the following,
tion comparing to the previous methods, it is not necessarily the mentioned challenges are discussed in more details.
the best PWL approximation.
In addition to the quadratic cost functions of generators, A. PWL Approximation for Convex Quadratic Functions
there are more instances of nonlinear functions in power In this section, the dimension of the problem is reduced to
system studies, many of which could be a multi-variable one for illustrative purposes. Assume that f (x) is a convex
function. The problem of fitting a PWL approximation to a quadratic function of the form
multi-dimensional set of data (possibly obtained by evaluating
a nonlinear function at different points) has been studied f (x) = ax2 + bx + c , a > 0, x ≤ x ≤ x (3)
before. The least-squares and Neural Network methods are
among the most-used approaches in curve-fitting. In order 1) Mathematical Background: It is obvious that the most
to find the best PWL approximation with a fixed number interesting PWL approximation of f is the one with the
of segments available, Magnani and Boyd have proposed a fewest linear parts and highest accuracy. It is also clear that
convex PWL fitting technique which is based on the so-called having more accuracy requires greater number of linear parts.
Max-Affine function [7]. This method is capable of providing However, if the PWL approximation is going to be used in,
PWL approximations for multivariate functions. This approach for example, an optimization problem, fewer segments is more
has specific application in convex optimization. interesting. Therefore, practical purposes limit the maximum
To the best of authors’ knowledge, this is the first time number of segments one can use in the PWL approximation.
in power system studies that the Max-Affine functions are Now, the question that remains is “having a limit on s, what
used for PWL applications. The capabilities of this method is the best fˆ(x)?”
are shown through numerical examples. This is the starting In order to answer the above question, it should be recalled
point for the numerous possibilities of these family of PWL that fˆ(x) can be expressed in more compact form as (referred
approximations in mathematical programming and optimiza- to as Max-Affine function in [7])
tion of power systems. The rest of the paper is organized as fˆ(x) = max {αi x + βi } (4)
follows. 1≤i≤s
In Section II, the PWL techniques are reviewed and illus- which surprisingly has no constraints on the subspaces on
trative examples are presented. In Section III the application which the linear approximations are defined. The best fˆ(x),
of the PWL technique in unit commitment problem and its with fixed s and m point-wise function evaluations, can then
impact on the solution quality and computational efficiency be obtained using the following least-squares problem:
are evaluated. The paper is concluded by highlighting the main
m 
contributions and findings of the present study. X 2
min max {αi xk + βi } − f (xk ) (5)
αi ,βi 1≤i≤s
II. P IECEWISE L INEAR A PPROXIMATION OF k=1
M ULTIVARIATE F UNCTIONS
Unfortunately, this problem is not convex [7]. However, an
Assume a function of multi variables, say f (x) with x ∈ efficient method is proposed in [7] to find the solution of (5).
Rn , is defined within a bound on x, say x ∈ D = {x|x ≤ This method is based on choosing the initial subspaces (i.e. Di )
x ≤ x}. Depending on the level of precision required, f (x) and updating them iteratively to find the best possible mesh
can be approximately expressed as PWL functions over small on D. Beside that algorithm, there are commercial solvers
sub-intervals inside D. For instance, considering s intervals, capable of handling these types of problems, which are usu-
one can derive the following as a PWL approximation of f (x): ally categorized as non-smooth problems. Some instances are
CONOPT, MINOS, LGO and IPOPT, all available in GAMS
 T
 a1 x + b1 , x ∈ D1
under the option “nonlinear programming with discontinuous

 aT2 x + b2 , x ∈ D2

fˆ(x) = .. (1) derivatives” [8]. The discussion on the algorithms for solving


 . non-smooth optimization problems is beyond the scope of this
 T
as x + bs , x ∈ Ds paper and the reader is referred to the software user manual.
in which ai ∈ Rn , bi ∈ R and the following holds:
[ \ 2) Numerical Example: Here, a numerical example for
Di = D and Di = ∅ (2) PWL approximation of a cost function for a thermal generation
1≤i≤s 1≤i≤s unit is presented. The coefficients in (3) are assumed to be
and on the borders of sequential Di , the linear segments are a = 0.9, b = 10, c = 200, x = 10, x = 200. Assume that
connected, which means that fˆ(x) is continuous. only two segments are allowed, i.e. s = 2. Figure 1 shows
Although (1) sounds interesting, there are real challenges the original quadratic function, the PWL upper approximation
in deriving fˆ(x). The first challenge is how to mesh D into (e.g. used in [3]) obtained by halving the space, and the
its subspaces Di . The second possible challenge is how to PWL max-affine approximation. The following remarks are
choose the smallest s while still maintaining a good accuracy observed:
 Table I
• Halving the available space is not the optimum way of AVERAGE OF THE R ELATIVE A BSOLUTE E RRORS FOR D IFFERENT VALUES
meshing (Points A and B in Fig. 1 are not equal). OF s O BTAINED BY T WO PWL TECHNIQUES
• The PWL upper approximation has no intersection with
the original function within the intervals but on the two Number of Partitions (s)
PWL Technique
1 2 3 4
ends.
Average of relative Upper Approximation 95.3 31.8 15.5 9.45
• The PWL max-affine approximation has 2 intersections absolute errors (%) Max-Affine Functions 90.5 20.1 9.2 1.3
with the original function within each interval (P1 , P2 in
the first interval and P3 , P4 in the second one in Fig. 1).
• The average of the relative absolute errors for the PWL
It is obvious that the PWL max-affine approximation is
upper approximation is 32% while this value for the PWL
more efficient than the other method. Moreover, the method
max-affine approximation is 20%.
is able to decide the optimal length of the intervals. Table
I shows the average of relative absolute errors for different
4
x 10
4
values of s. Also, Figs. 2 and 3 depict the histograms of
3.5
PWL max-affine the relative absolute errors for different values of s using the
Original quadratic
3 PWL upper approx. P
4
PWL max-affine and upper approximations, respectively. As
can be seen, the average relative error for the PWL max-affine
Cost ($/hour)

2.5

2
approximation is significantly less than the same values for the
1.5
A P PWL upper approximation.
3
1
B
0.5 P
0
2 III. A PPLICATIONS
P
1
-0.5
20 40 60 80 100 120 140 160 180 200 A. Min-Max Optimization
P (MW)
If the nonlinear function happens to be in the objective of
Figure 1. Comparison between the PWL upper and max-affine approximation a minimization programming, the PWL max-affine approxi-
techniques. mation leads to a linear reformulation of the objective. This
s=1
type of problems is referred to as Min-Max optimization.
s=2
0.8 0.6
Mathematically, it is described as
Normalized Density

0.6 0.45

0.4 0.3
min max {αiT x + βi } (6)
x∈D 1≤i≤s
0.2 0.15

0 0
By introducing a new variable, z = max{αiT x+βi }, the above
0.8
15 77 139 201 263 324 386 448
0.26
3 17 30 44 57 71 84 98
problem is reformulated as
s=3 s=4

min z (7a)
Normalized Density

0.6 0.20

x∈D
0.4 0.13

0.2 0.07 subject to z ≥ αiT x + βi , i = 1, . . . , s. (7b)

0 0.00
2 12 22 31 41 50 60 69 0.3 1.4 2.5 3.6 4.7 5.9 7.0 8.1 Therefore, by introducing one extra variable and s extra
Relative Error (%) Relative Error (%)
inequalities, the problem is reformulated as a linear pro-
Figure 2. Histograms of relative errors between the linearized and original gramming (assuming other constraints to be linear). This has
functions obtained using the PWL max-affine approximation. tremendous applications in power system optimization. As
an example, this method is applied to the unit commitment
s=1 s=2
0.15 0.3
problem in the following.
Normalized Density

0.12 0.24

0.09 0.18
B. Unit Commitment Problem
0.06 0.12

0.03 0.06 In this section, the linearization technique is applied to

0 0
4 18 32 46 60 75 89 103
the problem of unit commitment to show the impact of the
8 39 71 102 133 165 196 227

0.45 0.5
linearization on computational efficiency and quality. The unit
s=3 s=4
commitment problem is formulated here as follows (without
Normalizzed Density

0.36 0.4

0.27 0.3 loss of generality, some of the constraints are not considered
0.18 0.2 here for simplicity).
0.09 0.1

0 0
2 11 19 28 37 45 54 63 1 7 13 19 25 31 37 43
X X

Relative Error (%) Relative Error (%) Minimize CT = zi,h + SDni,h + SUpi,h (8)
i∈Ng h∈Nt
Figure 3. Histograms of relative errors between the linearized and original
functions obtained using the PWL upper approximation. subject to the following operational constraints:
 Table II
1) Active Power Flow Equations: U NITS ’ S CHEDULES O BTAINED BY MIQP FOR THE S IX -B US S YSTEM .
X
Pi,h − Pdi,h = B̄ij δj (9) Hour 1 2 3 4 5 6 7 8
j∈Nb G1 104.6 100 100 100 100 100 103.7 106
2) Line Flow Limits: G2 - - - - - - - -
X G6 84.2 78.1 71.1 66.8 67.2 73 83.2 85.5
−PiL ≤ D̄i,j δj ≤ PiL , i ∈ Nl (10) Hour 9 10 11 12 13 14 15 16
j∈Nb G1 110.9 123.1 146.5 125 128.3 129.1 131.9 135.6
G2 - - - 29.5 32.8 33.6 36.4 40.1
3) Generation Limits: G6 90.4 100 100 100 100 100 100 100
Hour 17 18 19 20 21 22 23 24
Pimin ui,h ≤ Pi,h ≤ Pimax ui,h (11)
G1 135.7 130.7 130.3 125.7 125.7 123.2 115.9 115.7
4) Shutdown/Startup Costs: G2 40.2 35.2 34.8 30.2 30.2 27.7 - -
G6 100 100 100 100 100 100 95.4 95.2
SUpi,h ≥ (ui,h − ui,h−1 )CiSUp , SUpi,h ≥ 0 (12)
Table III
SDni,h ≥ (ui,h−1 − ui,h )CiSDn , SDni,h ≥ 0 (13) U NITS ’ S CHEDULES O BTAINED BY MILP FOR THE S IX -B US S YSTEM .

5) Ramp Limits: Hour 1 2 3 4 5 6 7 8

G1 100 100 100 100 100 100 100 100
Pi,h − Pi,h−1 ≤ [ui,h − ui,h−1 ]PiSUp + ui,h−1 RiUp G2 - - - - - - - -
+ [1 − ui,h ]Pimax G6 88.9 78.1 71.1 66.8 67.2 73 86.9 91.5
Hour 9 10 11 12 13 14 15 16
(14)
G1 101.4 123.1 146.5 124 124 124 124 130
G2 - - - 30.5 37.1 38.6 44.3 46
Pi,h−1 − Pi,h ≤ [ui,h−1 − ui,h ]PiSDn + ui,h RiDn G6 100 100 100 100 100 100 100 100
Hour 17 18 19 20 21 22 23 24
+ [1 − ui,h−1 ]Pimax
G1 130 124 124 124 124 124 111.2 110.9
(15) G2 46 42 41.2 31.9 31.8 26.8 - -
G6 100 100 100 100 100 100 100 100
6) System Reserve:
X
ui,h Pimax − Pi,h ≥ PhRes


(16)
i∈Ng problem, for which there are efficient commercial solvers
7) Auxiliary Constraints: These constraints correspond to available, e.g. CPLEX [10]. All the problems are formulated
the PWL of the quadratic terms in the objective. in GAMS [8] and solved using CPLEX [10] in this paper.
For the PWL procedure, four segments are chosen, i.e.
zi,h ≥ αi,s Pi,h + βi,s ui,h (17) s = 4. Table II shows the units’ schedules obtained con-
Note that when ui,h = 0, it also follows that Pi,h = 0. sidering the exact quadratic cost functions (MIQP). Table III
Therefore, all the s inequalities turn to be zi,h ≥ 0, which the shows the units’ schedules for the six-bus system obtained
solver chooses the zero value. It is trivial to analyze the case using the linearized objective (MILP). As can be seen, the unit
of ui,h = 1. status (“on”/”off”) for both methods are identical. Also, the
differences between the committed generations are minor. The
C. Numerical Results computational efficiency of the two methods are compared in
In this section, the performance of the PWL approximation Table IV (obtained using an Intel Core i7-2600 CPU @ 3400
(i.e. the quality of the solution and the efficiency of the MHz). The “relative gap” is defined as the relative gap between
solution process) is presented through two examples. A Six- the best objective achieved up to the current iteration and the
bus system and the IEEE 118-bus test system are employed best lower bound. For the MIQP case, the solver could not
here for which the unit commitment problem is solved. The reach to the zero relative gap within the time limit of 1000
system data can be found in [9]. In the six-bus system, s, while for the MILP case, the proven optimal solution has
there are 3 generators and 7 branches. For the IEEE 118- been achieved within a few seconds. The objective value for
bus system, there are 54 generators and 186 branches. Two the six-bus system obtained using the MILP is higher than the
approaches are used to solve the unit commitment problem. one obtained by the MIQP. On the other hand, this is the other
In the first approach, the problem is formulated using the way around for the IEEE 118-bus system. These results reveal
original quadratic cost functions, which leads to a mixed- that both methods would come up with approximately same
integer quadratic programming (MIQP) problem. There are objective values.
commercial solvers capable of handling MIQP problems, It is worthwhile to compare the problem size with the
e.g. CPLEX [10]. In the second approach, the problem is method proposed in [3]. The number of extra variables re-
formulated using the PWL cost functions and associated extra quired for PWL approximation for each generator in [3] is
constraints, as given in Section III-B7. This leads to an MILP s + 1, while in the min-max formulation, only one extra
 Table IV
C OMPUTATIONAL E FFICIENCY OF MIQP AND MILP [4] A. Frangioni, C. Gentile, and F. Lacalandra, “Tighter approximated
MILP formulations for unit commitment problems,” IEEE Trans. Power
Syst., vol. 24, no. 1, pp. 105–113, Feb. 2009.
System Parameter MIQP MILP [5] A. Viana and J. P. Pedroso, “A new MILP-based approach for unit
CPU Time (s) 0.63 0.13 commitment in power production planning,” Int. Jour. Elect. Power &
Six-bus Objective Value ($/h) 153880 153975 Energy Syst., vol. 44, no. 1, pp. 997–1005, Jan. 2013.
Relative Gap 0 0 [6] L. Wu, “A tighter piecewise linear approximation of quadratic cost
curves for unit commitment problems,” IEEE Trans. Power Syst., vol. 26,
CPU Time (s) 1000 7.8
no. 4, pp. 2581–2583, Nov. 2011.
118-bus Objective Value ($/h) 650433 650387 [7] A. Magnani and S. Boyd, “Convex piecewise-linear fitting,” Optimiza-
Relative Gap 0.026 0 tion and Engineering, vol. 10, no. 1, pp. 1–17, 2009.
[8] “The General Algebraic Modeling System (GAMS).” [Online].
Available: http://www.gams.com
[9] Y. Fu, M. Shahidehpour, and Z. Li, “Security-constrained unit commit-
variable is needed. In addition, the number of constraints ment with AC constraints*,” IEEE Trans. Power Syst., vol. 20, no. 3,
for each generator in [3] is 2s + 2, while in the min-max pp. 1538–1550, Aug. 2005.
[10] “IBM ILOG CPLEX optimization studio.” [Online]. Avail-
formulation, only s constraints suffice. able: http://www-01.ibm.com/software/integration/optimization/cplex-
The applications of the PWL approximation is not limited to optimization-studio/
the unit commitment problem. There are many other problems [11] H. Ahmadi and H. Ghasemi, “Probabilistic optimal power flow incor-
porating wind power using point estimate methods,” in 10th Int. Conf.
in power system optimization which have a nonlinear objective Envir. Elec. Eng. (EEEIC), Rome, Italy, May 2011, pp. 1–5.
function subject to some linear constraints. For example, op- [12] M. Khanabadi, H. Ghasemi, and M. Doostizadeh, “Optimal transmission
timal power flow with objectives such as minimizing the cost switching considering voltage security and N-1 contingency analysis,”
IEEE Trans. Power Syst., vol. Early Access, no. 99, pp. 1–9, 2012.
or active losses is one good example, especially when there
is a need for multiple run [11]. As another instance, optimal
transmission line switching for congestion management could
be another application [12].

IV. C ONCLUSION
The nonlinear functions appearing in the objective function
of minimization problems are shown to be efficiently lineariz-
able using a piecewise linearization (PWL) technique. The
superiority of this method over existing PWL techniques are
demonstrated through examples. The main advantages of the
introduced approach can be summarized as follows:
• Higher accuracy in linearization is achieved by applying
the max-affine PWL technique.
• The size of the subspaces (on which the linear approxima-
tions are defined) are selected optimally by this method.
• The method is able to linearize a multivariate function in
multi-dimensional space.
• If the nonlinear function is in the objective, the advantage
can be taken of by minimizing a linear function subject
to a few affine inequalities.
• Significant saving in computation time can be achieved
when transforming a mixed-integer nonlinear program-
ming problem (with only the objective being nonlinear)
to a linear version.
As future work, further research is undertaken to apply this
competent technique to other areas of power system optimiza-
tion.

R EFERENCES
[1] Y.-C. Wu, A. S. Debs, and R. E. Marsten, “A direct nonlinear predictor-
corrector primal-dual interior point algorithm for optimal power flows,”
IEEE Trans. Power Syst., vol. 9, no. 2, pp. 876–883, May 1994.
[2] A. A. Sousa, G. L. Torres, and C. A. Cañizares, “Robust optimal power
flow solution using trust region and interior-point methods,” IEEE Trans.
Power Syst., vol. 26, no. 2, pp. 487–499, May 2011.
[3] M. Carrion and J. M. Arroyo, “A computationally efficient mixed-integer
linear formulation for the thermal unit commitment problem,” IEEE
Trans. Power Syst., vol. 21, no. 3, pp. 1371–1378, Aug. 2006.