Chapter 7

Optimal Dispatch of Generation Part I

7.1 Introduction

When generators are interconnected with various loads, the generation capacity is usually
much larger than the loads, hence the allocation of loads on the generators can be varied.
It is important to reduce the cost of electricity as much as possible, hence the "optimal"
division of load (meaning the least expensive) is desired. The most expensive cost of
generation is the fuel cost. Other expenses include labor, repairs and maintenance.
These are "economical" factors. Other factors (not considered in this presentation)
include such factors as security, stability and aesthetics.


7.2 Nonlinear Function Optimization

Unconstrained Parameter Optimization

The necessary condition for the function
( ) ( )
1 2
, , ,
n
f x f x x x K to have a minimum is
obtained by setting the derivatives 0, 1,2, ,
i
f
i n
x

K . Or using the operator,

0 f where
1 2
, , ,
n
f f f
f
x x x
_



,
K is known as the gradient of f . The term
H associated with second derivative:
2
i j
f
H
x x

results in a symmetric matrix called the
Hessian matrix of the function
( )
f x . Suppose
( )
0 f x where
( )
1 2
, , ,
n
x x x x K is a
local extremum. For this to be a minimum, the Hessian matrix should be positive
definite. This can be checked by finding the eigenvalues of the Hessian matrix, all of
which should be positive at x . In the process of doing this, two functions from Matlab
will be very useful, they are shown below. If x Ab , then the solution is found in
Matlab as \ x A b . Also, the eigenvectors of H are found using the Matlab function
eig(H).

Example 7.1

Find the minimum of
( )
2 2 2
1 2 3 1 2 2 3 1 2 3
2 3 8 16 32 110 f x x x x x x x x x x x + + + + + .
Finding the first derivatives and equate to zero results in three linear algebraic equations
thus:
2

1 2
1
1 2 3
2
2 3
3
2 8 0
4 16 0
6 32 0
f
x x
x
f
x x x
x
f
x x
x

+ +

or:

1
2
3
2 1 0 8
1 4 1 16
0 1 6 32
x
x
x
1 1 1
1 1 1

1 1 1
1 1 1
] ] ]

Using Matlab we have:

A=[2 1 0; 1 4 1; 0 1 6];
b=[8;16;32];
x=A\b
eigenvalues=eig(A)

x =
3.0000
2.0000
5.0000
eigenvalues =
1.5505
4.0000
6.4495

Thus the eigenvalues are all positive, thus x=(3, 2, 5) is a minimum (in this case a global
minimum since there is only one extremum).


7.2.1 Constrained Parameter Optimization:
Equality Constraints

Here we minimize
( )
f x subject to the constraints
( )
0, 1,2, ,
i
g x i k K . Such
problems can be solved using Lagrange multipliers
i
thus we let
1
k
i i
i
f g

L . The
necessary conditions for constrained local minima of L are the following:

1
0
0
k
i
i
i i i i
i
i
f g
x x x
g

L
L

Note that the last equation is the same as the original constraints.

Example 7.2
3

Use the Lagrange multiplier method for solving constrained parameter optimizations to
determine the minimum distance from the origin of the xy-plane to a circle described by
( ) ( )
2 2
8 6 25 x y + . The distance to be minimized is given by
( )
2 2
, f x y x y + .

First we use Matlab to show a graphical display of the two curves from which it is clear
that the answer is at (4, 3). Note there is also a maximum at (12, 9).

wt=0:.01:2*pi;
z =8+j*6+ 5*(cos(wt) + j*sin(wt));
x=0:.01:12; y=6/8*x;
plot(real(z),imag(z), x, y), grid
xlabel('x'), ylabel('y')
axis([0 , 14, 0, 14])

0 2 4 6 8 10 12 14
0
2
4
6
8
10
12
14
x
y


Now we want to minimize
( )
, f x y subject to the constraints described by the circle
equation. We first form the Lagrange function thus:

( ) ( )
2 2
2 2
8 6 25 x y x y
1
+ + +
]
L
The necessary conditions for an extremum are found as
4

( ) ( )
( ) ( )
( ) ( )
2 2
2 2 16 0 or 2 1 16
2 2 12 0 2 1 12
8 6 25 0
x x x
x
y y or y
y
x y

+ +

+ +

L
L
L

Dividing the first two equations yields
3
4
y x . Using this in the third equation yields

2
25
25 75 0
16
x x +
Solving the quadratic gives 4 x and 12 x . Thus the extrema are at the points (4, 3)
with 1 , and (12, 9) with 3 . From the figure it is clear that the first point is the
minimum.

In many situations the equations are not only nonlinear but also cannot be solved in an
analytical manner, thus the solution has to be obtained by iteration. The most common
iteration method is the Newton-Raphson method. This is done below for the above
problem. Solve the first two equations for x and y in terms of :

8
1
6
1
x
y

+

Using these values in the third equation yields:

( ) ( ) ( )
( )
2
2
100 200
, 75 0
1
1
f f x y

+ 1
]
+
+

This equation can be solved by iteration as follows using the Newton-Raphson method.
Using the first order term in a Taylor series expansion we have:

( )
( )
( )
( )
( ) ( ) ( ) 1
and
k
k k k k
k
f
df
d

+
_

,

where
( )
( ) k
f is the residual at
( ) k
given by
( )
( )
( ) ( )
( )
,
k k k
f f x y . If this residual is
zero, the solution is exact, otherwise we would like the residual to satisfy a very small
convergence value. The iteration is started with an estimated value of and continued
till a specified accuracy. Once is known, x and y can be computed from the equations
above. One disadvantage to this method is the need to find the derivative of the function
which in this case is found to be:

( )
( ) ( ) ( )
3 2 3
200 200 200
1 1 1
df
d

+ + +

5
Having a value for
(0)
using the above equation we have
(0)
f

,
. Now we can
compute
(0)
x and
(0)
y hence we have a value for use in:
( ) ( ) ( )
2 2
(0)
( 0) (0)
8 6 25 f x y +
and now we can find
( )
( )
( )
0
(0)
0
f
df
d

_

,
. Thus we have a new value for
(1)
. The
process is repeated till the the error in evaluating f is less than some prescribed value.

The following Matlab program will perform the iteration:

clear;
iter = 0; % Iteration counter
Df = 10; % Error in Df is set to a high value
Lambda = 0.4; % Initial estimated value of Lambda
disp([' Iter Df J DLambda Lambda x
y'])
while abs(Df) >= 0.0001 % Test for convergence
iter = iter + 1; % No. of iterations
x = 8*Lambda/(Lambda + 1);
y = 6*Lambda/(Lambda + 1);
Df = (x- 8)^2 + (y - 6)^2 - 25; % Residual
J = 200*Lambda/(Lambda + 1)^3 - 200/(Lambda + 1)^2; % Derivative
Delambda =-Df/J; % Change in variable
disp([iter, Df, J, Delambda, Lambda, x, y])
Lambda = Lambda + Delambda; % Successive solution
end

Iter Df J DLambda Lambda x y
1.0000 26.0204 -72.8863 0.3570 0.4000 2.2857 1.7143
2.0000 7.3934 -36.8735 0.2005 0.7570 3.4468 2.5851
3.0000 1.0972 -26.6637 0.0411 0.9575 3.9132 2.9349
4.0000 0.0337 -25.0505 0.0013 0.9987 3.9973 2.9980
5.0000 0.0000 -25.0001 0.0000 1.0000 4.0000 3.0000

We repeat the program with an initial value of 2 :

clear;
iter = 0; % Iteration counter
Df = 10; % Error in Df is set to a high value
Lambda = -2; % Initial estimated value of Lambda
disp([' Iter Df J DLambda Lambda x
y'])
while abs(Df) >= 0.0001 % Test for convergence
iter = iter + 1; % No. of iterations
x = 8*Lambda/(Lambda + 1);
y = 6*Lambda/(Lambda + 1);
Df = (x- 8)^2 + (y - 6)^2 - 25; % Residual
J = 200*Lambda/(Lambda + 1)^3 - 200/(Lambda + 1)^2; % Derivative
Delambda =-Df/J; % Change in variable
6
disp([iter, Df, J, Delambda, Lambda, x, y])
Lambda = Lambda + Delambda; % Successive solution
end

Iter Df J DLambda Lambda x y
1.0000 75.0000 200.0000 -0.3750 -2.0000 16.0000 12.0000
2.0000 27.8926 76.9346 -0.3625 -2.3750 13.8182 10.3636
3.0000 8.1227 38.1258 -0.2131 -2.7375 12.6042 9.4531
4.0000 1.2823 26.9480 -0.0476 -2.9506 12.1013 9.0760
5.0000 0.0454 25.0682 -0.0018 -2.9982 12.0036 9.0027
6.0000 0.0001 25.0001 -0.0000 -3.0000 12.0000 9.0000

Now the solution converges to the "maximum" distance at (12, 9).


7.2.2 Constraint Parameter Optimization:
Inequality Constraints

In practice both equality and inequality constraints are needed. The problem is to
minimize the cost function

( )
1 2
, , ,
n
f x x x K
subject to the equality constraints

( )
1 2
, , , 0 1,2, ,
i n
g x x x i k K K
and inequality constraints

( )
1 2
, , , 0 1,2, ,
j n
u x x x j m K K
The Lagrange multipliers are extended to include the inequality constraints using the
multipliers as shown below:

1 1
k m
i i j j
i j
f g u

+ +

L
The resulting necessary conditions for constrained local minima of L are the following:

0 1, ,
0 1, ,
0 1, ,
0 & 0 1, ,
i
i
i
j
j
j j j
i n
x
g i k
u i m
u i m

>
K
K
K
K
L
L
L

Note that the second equation above is the equality constraints and the last two equations
pertain to the inequality constraints. Suppose
( )
1 2
, , ,
n
x x x K is a relative minimum. The
inequality constraint (third equation above) is said to be inactive if strict inequality holds
at
( )
1 2
, , ,
n
x x x K and 0
j
. On the other hand, when strict equality holds, the
7
constraint is active at this point (i.e. if the constraint
( )
1 2
, , , 0
j j n
u x x x K and
0
j
> ). This is known as the Kuhn-Tucker necessary condition.

Example 7.3

Solve example 7.2 with an additional inequality constraint defined below. The problem
is to find the minimum of the function
( )
2 2
, f x y x y + subject to one equality
constraint
( ) ( ) ( )
2 2
, 8 6 25 0 g x y x y + and one inequality constraint
( )
, 2 12 u x y x y + .

The unconstrained cost function is L as shown below:
( ) ( ) ( )
2 2
2 2
8 6 25 2 12 x y x y x y
1
+ + + + +
]
L
The resulting necessary conditions for constrained local minima of L are

( )
( )
( ) ( )
2 2
2 2 8 2 0
2 2 6 0
8 6 25 0
2 12 0
x x
x
y y
y
x y
x y

+ +

+ +

L
L
L
L

Eliminating from the first two equations results in

( )( )
2 4 1 8 0 x y + +
From the fourth condition we have
12 2 y x
Now substituting for y in the previous equation yields

4 4.8
1
x

+

And substituting for x in the fourth condition we get

4 2.4
1
y

+

Substituting for x and y in the third condition (equality constraint) results in an equation
in terms of :

2 2
4 4.8 4 2.4
8 6 25 0
1 1


+ +
_ _
+

+ +
, ,

from which we have the following equation

2
2 0.36 0 + +
Roots of this equation are at 0.2 and 1.8 . Substituting for these values of
in the expressions for x and y, the corresponding extrema are
8

( ) ( )
( ) ( )
, 5,2 for 0.2, 5.6
, 3,6 for 1.8, 12
x y
x y






Using Matlab these are graphically depicted below. Note that the solution is valid above
the line shown in the figure below, thus only two points satisfy the solution for extrema,
one is a minimum, the other a maximum. These are the points where the circle and the
line intersect, which are exactly the two points found above. The graphical
demonstration below further illustrates this example.

%Graphical Demonstration of Example 7.3

wt=0:.01:2*pi;
z =8+j*6+ 5*(cos(wt) + j*sin(wt));
x=0:.01:12; y=6/8*x;
y2=12-2*x;
plot(real(z),imag(z), x, y, x, y2), grid
xlabel('x'), ylabel('y')
axis([0 , 14, 0, 14])


a =[1 2 .36];
Lambda = roots(a)
X=(4*Lambda+4.8)./(1+Lambda)
Y=(4*Lambda+2.4)./(1+Lambda)
F=sqrt(X.^2+Y.^2)
Mindist=min(F)

Lambda =
-1.8000
-0.2000
X =
3.0000
5.0000
Y =
6.0000
2.0000
F =
6.7082
5.3852
Mindist =
5.3852
9
0 2 4 6 8 10 12 14
0
2
4
6
8
10
12
14
x
y