Ch- 2 : Economic Operation of Power System & Unit Commitment

exa 2.2 page 39

In [9]:
from __future__ import division
from sympy.mpmath import quad
#For equal incremental cost
L1=125 #MW
L2=100 #MW
#For equal sharing
L=(L1+L2)/2 #MW
#Change in cost Unit 1
f1 = lambda P1:0.2*P1+30 
dC1=quad(f1,[L1,L]) #Rs./hour
#Change in cost Unit 2
f2 = lambda P2 : 0.15*P2+40
dC2=quad(f2,[L2,L]) #Rs./hour
dCyearly=(dC1+dC2)*24*365 #Rs./year
print "Saving per year in economic load allocation = %0.2f Rs./year "%(dCyearly) 
#Answer in the textbook is not accurate.

Saving per year in economic load allocation = 239531.25 Rs./year

exa 2.3 page 43

In [10]:
L=400 #/MW#/total load
delPD=50 #MW#increase in demand
#dC1/dP1=0.2*P1+30
#dC2/dP2=0.15*P2+40
twoC1=0.2 ##from above equation
twoC2=0.15 ##from above equation
delP1_by_delPD=(1/twoC1)/(1/twoC1+1/twoC2) 
delP2_by_delPD=(1/twoC2)/(1/twoC1+1/twoC2) 
delP1=delP1_by_delPD*delPD #MW
print "Increase in generation of unit1 = %0.2f MW  "%delP1 
delP2=delP2_by_delPD*delPD #MW
print "Increase in generation of unit2 = %0.2f MW  "%delP2 
P1=L/2+delP1 #load on unit 1
print "Total load on unit1 = %0.2f MW "%P1 
P2=L/2+delP2 #load on unit 2
print "Total load on unit2 = %0.2f MW" %P2
print "Checking incremental cost :" 
dC1_by_dP1=0.2*P1+30 #Rs./MWh
print "Incremental cost of unit 1 = %0.2f Rs./MWh "%dC1_by_dP1 
dC2_by_dP2=0.2*P2+30 #Rs./MWh
print "Incremental cost of unit 2 = %0.2f Rs./MWh "%dC2_by_dP2 
print "Conclusion : Cost are same(Approximately)." 
#Note : Values calculated in the book are slightly wrong because of accuracy in calculation as compared to scilab accuracy.

Increase in generation of unit1 = 21.43 MW  
Increase in generation of unit2 = 28.57 MW  
Total load on unit1 = 221.43 MW 
Total load on unit2 = 228.57 MW
Checking incremental cost :
Incremental cost of unit 1 = 74.29 Rs./MWh 
Incremental cost of unit 2 = 75.71 Rs./MWh 
Conclusion : Cost are same(Approximately).

exa 2.5 page 47

In [11]:
from math import atan, cos
I1=0.8 #p.u.
I2=1 #p.u.
Za=0.04+1J*0.12 #p.u.
Zb=0.03+1J*0.1 #p.u.
Zc=0.03+1J*0.12 #p.u.
V=1 #p.u.
#Solution : 
V1=V+(I1+I2)*Za+I1*(Zb) #p.u.
V2=V+(I1+I2)*Za+I2*(Zc) #p.u.
P1=(I1*V1).real #p.u.
P2=(I2*V2).real #p.u.
fi1=atan((V1.imag)/(V1.real)) 
fi2=atan((V2.imag)/(V2.real)) 
print "Loss Coefficients are : "
B11=((Za.real)+(Zb.real))/(abs(V1)**2*cos(fi1)**2) #p.u.
print "B11 = %0.5f"%B11,"p.u." 
B22=((Za.real)+(Zc.real))/(abs(V2)**2*cos(fi2)**2) #p.u.
print "B22 = %0.4f"%B22,"p.u. " 
B12=((Za.real))/(abs(V1)*abs(V2)*cos(fi1)*cos(fi2)) #p.u.
print "B12 = %0.4f"%B12,"p.u." 
PL=P1**2*B11+P2**2*B22+2*P1*P2*B12 #p.u.
print "Transmission Loss = %0.6f p.u." %PL
#Note : Values calculated in the book are slightly wrong because of accuracy in calculation as compared to scilab accuracy.

Loss Coefficients are : 
B11 = 0.05827 p.u.
B22 = 0.0576 p.u. 
B12 = 0.0331 p.u.
Transmission Loss = 0.178800 p.u.

exa 2.7 page 53

In [12]:
from math import atan, degrees, cos, pi
Za=0.03+1J*0.09 #p.u.
Ia=1.5-1J*0.4 #p.u.
Zb=0.10+1J*0.30 #p.u.
Ib=0.5-1J*0.2 #p.u.
Zc=0.03+1J*0.09 #p.u.
Ic=1-1J*0.1 #p.u.
Zd=0.04+1J*0.12 #p.u.
Id=1-1J*0.2 #p.u.
Ze=0.04+1J*0.12 #p.u.
Ie=1.5-1J*0.3 #p.u.
V=1 #p.u.
base=100 #MVA
#Solution
#Currents of load
IL1=0.4 #p.u.
IL2=0.6 #p.u.
#Current distribution factors :
Na1=1; Na2=0 
Nb1=0.6; Nb2=-0.4 
Nc1=0 ;Nc2=1 
Nd1=0.4; Nd2=0.4 
Ne1=0.6 ;Ne2=0.6 
#Bus Voltages
V1=V+Ia*Za #p.u.
V2=V-Ib*Zb+Ic*Zc #p.u.
#Phase Angles
theta1=degrees(atan((Ia.imag)/(Ia.real)) )#degree
theta2=degrees(atan((Ic.imag)/(Ic.real)) )#degree
#Power Factors : 
cos_fi1=cos(atan((V1.imag)/(V1.real))-theta1*pi/180) #source 1 power factor
cos_fi2=cos(atan((V2.imag)/(V2.real))-theta2*pi/180) #source 2 power factor
print "Loss formula Coefficients in p.u. :"
B11=(Na1**2*(Za.real)+Nb1**2*(Zb.real)+Nc1**2*(Zc.real)+Nd1**2*(Zd.real)+Ne1**2*(Ze.real))/(abs(V1)**2*cos_fi1) #p.u.
print "B11 = %0.5f p.u "%B11 
B22=(Na2**2*(Za.real)+Nb2**2*(Zb.real)+Nc2**2*(Zc.real)+Nd2**2*(Zd.real)+Ne2**2*(Ze.real))/(abs(V2)**2*cos_fi2) #p.u.
print "B22 = %0.5f p.u" %B22
B12=(Na1*Na2*(Za.real)+Nb1*Nb2*(Zb.real)+Nc1*Nc2*(Zc.real)+Nd1*Nd2*(Zd.real)+Ne1*Ne2*(Zc).real)/(abs(V1)*abs(V2)*cos_fi1*cos_fi2*(cos(theta1*pi/180-theta2*pi/180))) #p.u.
print "B12 = %0.5f"%B12,"p.u " 
#Converting p.u. to actual value
print "Loss formula Coefficients in MW**-1 :"
B11=B11/base #MW**-1
print "B11 = %0.5f"%B11,"MW**-1" 
B22=B22/base #MW**-1
print "B22 = %0.5f"%B22,"MW**-1" 
B12=B12/base #MW**-1
print "B12 = %0.5f"%B12,"MW**-1" 
#Note : Values calculated in the book are slightly wrong because of accuracy in calculation as compared to scilab accuracy.

Loss formula Coefficients in p.u. :
B11 = 0.07877 p.u 
B22 = 0.07735 p.u
B12 = -0.00732 p.u 
Loss formula Coefficients in MW**-1 :
B11 = 0.00079 MW**-1
B22 = 0.00077 MW**-1
B12 = -0.00007 MW**-1

exa 2.8 page 54

In [13]:
#dC1/dP1=0.2*P1+22 #Rs./MWh
#dC2/dP2=0.15*P2+30 #Rs./MWh
B22=0; B12=0 #Because Loss is independent wrt P2
P1=100 #MW
PL=15 #MW
B11=PL/P1**2 #MW**-1
L1=1/(1-0.003*P1)#Penalty Factor plant 1
L2=1 #Penalty Factor of plant 2
lamda=60 
#lamda=dC1/dP1*L1=dC2/dP2*L2
#dC1/dP1*L1=dC2/dP2*L2
P2=((0.2*P1+22)*L1-30)/0.15 #MW
P=P1+P2-B11*P1**2 #MW#Total Load
print "Required generation at plant1 = %0.2f MW "%P1 
print "Required generation at plant2 = %0.2f MW" %P2
print "Total Load = %0.2f MW "%P 

Required generation at plant1 = 100.00 MW 
Required generation at plant2 = 200.00 MW
Total Load = 285.00 MW

exa 2.9 page 55

In [14]:
from sympy import symbols, solve
#dC1/dP1=0.2*P1+22 #Rs./MWh
#dC2/dP2=0.15*P2+30 #Rs./MWh
B22=0 ;B12=0 #Because Loss is independent wrt P2
P1=100 #MW
PL=15 #MW
B11=PL/P1**2 #MW**-1
L1=1/(1-0.003*P1) #Penalty Factor plant 1
L2=1 #Penalty Factor of plant 2
lamda=60 
#lamda=dC1/dP1*L1=dC2/dP2*L2
#dC1/dP1*L1=dC2/dP2*L2
P2=((0.2*P1+22)*L1-30)/0.15 #MW
P=P1+P2-B11*P1**2 #MW#Total Load
#dC1/dP1=dC2/dP2  neglecting transmission loss
#clear('P2') #for recalculation
#0.2*P1-0.15*P2-8=0 #/eqn(1)
#P1=0.75*P2+40 #P1+P2-B11*P1**2-P=0 #eqn(2)
#1.75*P2-B11*P1**2=P-40
#Eqn=[-B11 1.75 40-P] 
P22 = symbols('P22')
Eqn=-B11*P22**2+1.75*P22+40-P 
P22=solve(Eqn, P22)
P2=P22[0]#MW#neglecting higher value
P1=0.75*P2+40 #MW
#from sympy.mpmath import quad
dC1=quad(lambda P: 0.2*P+22,[100,P1]) #Rs.#Additional Cost plant1
dC2=quad(lambda P: 0.15*P+30,[200,P2]) #Rs.#Decreased Cost  plant2
dC=dC1+dC2 #Rs./hour#Net change in cost
print "Taking transmission loss in account, Net saving per hour in fuel cost = %0.2f Rs./hour "%dC 
#Note : Values calculated in the book are slightly wrong because of accuracy in calculation as compared to scilab accuracy.

Taking transmission loss in account, Net saving per hour in fuel cost = 854.74 Rs./hour

exa 2.10 page 55

In [15]:
B11=0.001 #MW**-1
B22=0.0024 #MW**-1
B12=-0.0005 #MW**-1
#dC1/dP1=0.8*P1+16 #Rs./MWh
#dC2/dP2=0.08*P2+12 #Rs./MWh
lamda=20 
#Iterations for calculating value
P1 = range(0,10)
P1[0]=0 
P2 = range(0,10)
P2[0]=0 
for i in range(1,10):
    P1[i] =( 0.2+0.001*P2[i-1])/0.006 
    P2[i] = (0.4+0.001*P1[i])/0.0088 
    if P1[i]==P1[i-1]:
        break 
    

P1=P1[i] #MW
print "Generation P1 = %0.2f MW" %P1
P2=P2[i] #MW
print "Generation P2 = %02.f MW "%P2 
PL=B11*P1**2+2*B12*P1*P2+B22*P2**2 #MW
print "Transmission Loss = %0.2f MW "%PL 
Pr=P1+P2-PL #MW
print "Received Power = %0.2f MW "%Pr 

Generation P1 = 41.70 MW
Generation P2 = 50 MW 
Transmission Loss = 5.69 MW 
Received Power = 86.20 MW

exa 2.11 page 56

In [16]:
from sympy import symbols
#C1=561+7.92*P1+0.001562*P1**2 #Rs./hour
#C2=310+7.85*P2+0.00194*P2**2 #Rs./hour
a1=561 ;a2=310 
b1=7.92 ;b2=7.85 
c1=0.001562 ;c2=0.00194 
ce=c1*c2/(c1+c2) 
be=ce*(b1/c1+b2/c2) 
ae=a1-b1**2/4/c1+a2-b2**2/4/c2+be**2/4/ce 
print "Coefficients are : " 
print "ae = ",(ae)," & be = ",(be)
print "ce = ",ce
PT = symbols('PT') 
CT = round(ae,2)+round(be,2)*PT+round(ce,4)*PT**2
print "Cost Characteristics : ",CT
#print "CT=870.753+7.8888*PT+0.0008653*PT**2" 

Coefficients are : 
ae =  870.650199886  & be =  7.88877784123
ce =  0.000865299828669
Cost Characteristics :  0.0009*PT**2 + 7.89*PT + 870.65

exa 2.12 page 57

In [17]:
from sympy import symbols,solve
#C1=7700+52.8*P1+5.5*10**-3*P1**2 #Rs./hour
#C2=2500+15*P2+0.05*P2**2 #Rs./hour
a1=7700 ;a2=2500 
b1=52.8 ;b2=15 
c1=5.5*10**-3 ;c2=0.05 
P1, P2 = symbols('P1 P2')
dC1bydP1=52.8+2*5.5*10**-3*P1 
dC2bydP2=15+2*0.05*P2 
print "For 1200 MW Load :" 
P=1200 #MW
#Let loads of unit are P1 & 1200-P1
#Economical Loading dC1/dP1=dC2/dP2
eqn=52.8+2*5.5*10**-3*P1-15-2*0.05*(1200-P1) 
P1=solve(eqn, P1)[0] #MW
P2=P-P1 #MW
print "P1 = %0.2f MW "%P1 
print "P2 = %0.2f MW "%P2 
print "For 900 MW Load :" 
P=900 #MW
#clear('P1','P2') 
P1, P2 = symbols('P1 P2')
#Let loads of unit are P1 & 900-P1
#Economical Loading dC1/dP1=dC2/dP2
eqn=52.8+2*5.5*10**-3*P1-15-2*0.05*(900-P1) 
P1=solve(eqn, P1)[0] #MW
P2=P-P1 #MW
print "P1 = %0.2f MW "%P1 
print "P2 = %0.2f MW "%P2 
print "For 500 MW Load :" 
P=500 #MW
#clear('P1','P2') 
P1, P2 = symbols('P1 P2')
#Let loads of unit are P1 & 500-P1
#Economical Loading dC1/dP1=dC2/dP2
eqn=52.8+2*5.5*10**-3*P1-15-2*0.05*(500-P1) 
P1=solve(eqn, P1)[0] #MW
P2=P-P1 #MW
#Minimum load is 200MW
if P1<200:
    P2=P1+P2
    P1=0 
print "P1 = %0.2f MW "%P1 
print "P2 = %0.2f MW "%P2 
C=(2500+15*P2+0.05*P2**2)*10 #Rs.#Operating cost for 10 hour
print "Operating cost for 10 hour = %0.2f Rs. "%C 
print "Other option : " 
P1=200 #MW
P2=300 #MW
print "P1 = %0.2f MW "%P1 
print "P2 = %0.2f MW "%P2 
C1=7700+52.8*P1+5.5*10**-3*P1**2 #Rs./hour
C2=2500+15*P2+0.05*P2**2 #Rs./hour
C=10*(C1+C2) #Rs.#Operating cost for 10 hour
print "Operating cost for 10 hour = %0.2f Rs. "%C 

For 1200 MW Load :
P1 = 740.54 MW 
P2 = 459.46 MW 
For 900 MW Load :
P1 = 470.27 MW 
P2 = 429.73 MW 
For 500 MW Load :
P1 = 0.00 MW 
P2 = 500.00 MW 
Operating cost for 10 hour = 225000.00 Rs. 
Other option : 
P1 = 200.00 MW 
P2 = 300.00 MW 
Operating cost for 10 hour = 299800.00 Rs.

exa 2.13 page 58

In [18]:
#C1=2000+20*P1+0.05*P1**2 #Rs./hour
#C2=2750+26*P2+0.03091*P2**2 #Rs./hour
P1=350 #MW
P2=550 #MW
C1=2000+20*P1+0.05*P1**2 #Rs./hour
C2=2750+26*P2+0.03091*P2**2 #Rs./hour
C=C1+C2 #Rs./hour
print "(a) Total Cost = %0.2f Rs./hour "%C 
P=P1+P2 #MW#Total Load
P1, P2 = symbols('P1 P2')
dC1bydP1=20+2*0.05*P1 
dC2bydP2=26+2*0.03091*P2 
print "(b) For Economic Scheduling"
#dC1/dP1=dC2/dP2 for economic sheduling
#Let loads of unit are P1 & P-P1
eqn=20+2*0.05*P1-26-2*0.03091*(P-P1) 
P1=solve(eqn, P1)[0] #MW
P2=P-P1 #MW
print "Loads P1 & P2 = %0.2f & %0.2f MW"%(P1, P2)
C1=2000+20*P1+0.05*P1**2 #Rs./hour
C2=2750+26*P2+0.03091*P2**2 #Rs./hour
Cnew=C1+C2 #Rs./hour
print "Total Cost = %0.2f Rs./hour "%Cnew 
saving=C-Cnew #Rs./hour
print "Total saving  %0.2f Rs./hour "%saving 
Lt=P1-350 #MW#Tie line load
print "Tie line load from Plant1 to Plant2 = %0.2f MW  "%Lt

(a) Total Cost = 41525.28 Rs./hour 
(b) For Economic Scheduling
Loads P1 & P2 = 380.90 & 519.10 MW
Total Cost = 41448.00 Rs./hour 
Total saving  77.28 Rs./hour 
Tie line load from Plant1 to Plant2 = 30.90 MW

exa 2.14 page 59

In [19]:
from sympy import symbols, solve
#C=5000+450*P+0.5*P**2 #Rs./hour
e1=+2 #%#error
e2=-2 #%#error
P=200 #MW#Total Load
#Considering error
P1, P2 = symbols('P1 P2')
C1=(5000+450*P+0.5*P1**2)*0.98 #Rs./hour
C2=(5000+450*P+0.5*P2**2)*1.02 #Rs./hour
#Let loads of unit are P1 & P-P1
#dC1/dP1=dC2/dP2 for economic sheduling
eqn=450*0.98+2*0.5*P1*0.98-450*1.02-2*0.5*(P-P1)*1.02 
P1=solve(eqn, P1)[0] #MW
P2=P-P1 #MW
#if no instrumention error
C1=(5000+450*P1+0.5*P1**2)*0.98 #Rs./hour
C2=(5000+450*P2+0.5*P2**2)*1.02 #Rs./hour
C=C1+C2 #Rs./hour
#Due to intrumentation error
P1=P/2 #MW
P2=P/2 #MW
C1=(5000+450*P1+0.5*P1**2)*0.98 #Rs./hour
C2=(5000+450*P2+0.5*P2**2)*1.02 #Rs./hour
Cerr=C1+C2 #Rs./hour
Cextra=Cerr-C #Rs,/hour
print "Extra operating cost = %0.2f Rs./hour "%Cextra 

Extra operating cost = 121.00 Rs./hour

exa 2.15 page 59

In [20]:
from sympy import symbols, solve
import numpy as np
P1, P2, P3 = symbols('P1 P2 P3')
Q1=0.002*P1**2+0.86*P1+20 #tons/hour
Q2=0.004*P2**2+1.08*P2+20 #tons/hour
Q3=0.0028*P3**2+0.64*P3+36 #tons/hour
Pmax=120 #MW
Pmin=36 #MW
P=200 #MW
C=500 #Rs./ton
#C1=C*Q1 C2=C*Q2 C3=C*Q3 #Rs./ton
dC1bydP1=2*P1+430 #Rs./hour
dC2bydP2=4*P2+540 #Rs./hour
dC3bydP3=2.8*P3+320 #Rs./hour
#P1+P2+P3=P
A1=[1 ,1 ,1] #Coefficient Matrix
B1=[P] #Coefficient Matrix
#For minimal cost above 3 equation should be equal
#eqn1=2*P1-4*P2+430-540 
#eqn2=4*P2-2.8*P3-320+540 
A2=[0 ,4 ,-2.8] #Coefficient Matrix
B2=[-540+320] #Coefficient Matrix
#eqn3=-2*P1+2.8*P3+320-430 
A3=[-2 ,0 ,2.8] #Coefficient Matrix
B3=[430-320] #Coefficient Matrix
#solving by matrix method
A=np.mat([A1[:] ,A2[:], A3[:]]) #Coefficient Matrix
B=[B1 ,B2 ,B3] #Coefficient Matrix
X=A**-1*B #Solution Matrix
P1=X[0] #MW
P2=X[1] #MW
P3=X[2] #MW
Pmax=120 #MW
Pmin=36 #MW
if P2<Pmin:
   P2=Pmin #MW     

#P1+P3=P-P2#eqn(4)
A1=[1 ,1] #Coefficient Matrix
B1=[P-P2] #Coefficient Matrix
#eqn3=-2*P1+2.8*P3+320-430 
A2=[-2, 2.8] #Coefficient Matrix
B2=[430-320] #Coefficient Matrix
#solving by matrix method
A=np.mat([A1[:], A2[:]]) #Coefficient Matrix
B=np.mat([B1[:], B2[:]]) #Coefficient Matrix
X=A**-1*B #Solution Matrix
P1=X[0] #MW
P3=X[1] #MW
print "According to optimum scheduling, Load distriution is :" 
print "P1 = %0.2f MW "%P1 
print "P2 = %0.2f MW "%P2 
print "P3 = %0.2f MW "%P3 

According to optimum scheduling, Load distriution is :
P1 = 72.75 MW 
P2 = 36.00 MW 
P3 = 91.25 MW

exa 2.16 page 60

In [21]:
L=30 #MW
#I=(32+32*L+1.68*L**2)*10**5 
t1=18 #/hours
t2=6 #/hours
#Full load 18 hours
I1=(32+32*L+1.68*L**2)*10**5*t1 #kJ
#Half load 6 hours
I2=(32+32*L/2+1.68*(L/2)**2)*10**5*t2
I=I1+I2 #kJ
print "(a) Heat input per day = %0.2e kJ "%I 
E=L*t1+L/2*t2 #MWh#/Energy produced in 24 hours
Lu=E/(t1+t2) #MW
Inew=(32+32*Lu+1.68*Lu**2)*10**5*(t1+t2) #kJ
saving=I-Inew #/kJ
saving=saving/(E*1000) #kJ/kWh
print "(b) Saving in heat per kWh of energy = %0.2f kJ/kWh " %saving

(a) Heat input per day = 5.04e+09 kJ 
(b) Saving in heat per kWh of energy = 270.00 kJ/kWh

exa 2.17 page 61

In [22]:
import numpy as np
P=800 #MW(Total Load)
#Using Variable for Cost Curve Equation
P1, P2, P3 = symbols('P1 P2 P3')
#Cost Curve Equation
C1=450+6.5*P1+0.0013*P1**2 #Rs./hour
C2=300+7.8*P2+0.0019*P2**2 #Rs./hour
C3=80+8.1*P3+0.005*P3**2 #Rs./hour
#Part(a) is not computational
#Part (b)
dC1BydP1=6.5+2*0.0013*P1 #Rs./MWh#/eqn(1)
dC2BydP2=7.8+2*0.0019*P2 #Rs./MWh#/eqn(2)
dC3BydP3=8.1+2*0.005*P3 #Rs./MWh#/eqn(3)
#P1+P2+P3=P #MW#/eqn(4)
A1=[1 ,1 ,1] #Coefficient Matrix
B1=[800] #Coefficient Matrix
#Equating eqn(1) & (2)
A2=[2*0.0013, -2*0.0019, 0] #Coefficient Matrix
B2=[7.8-6.5] #Coefficient Matrix
#Equating eqn(2) & (3)
A3=[0 ,2*0.0019, -2*0.005] #Coefficient Matrix
B3=[8.1-7.8] #Coefficient Matrix
#Solution By Matrix method
A=np.mat([A1[:], A2[:], A3]) #Coefficient Matrix
B=np.mat([B1[:], B2[:], B3[:]]) #Coefficient Matrix
X=A**-1*B #Solution Matrix
P1=X[0] #MW
P2=X[1] #MW
P3=X[2] #MW
print "(b) According to optimum scheduling, Load distriution is :" 
print "P1 = %0.2f MW" %P1
print "P2 = %0.2f MW" %P2
print "P3 = %0.2f MW" %P3
#Part(c)
print "(c) Optimum scheduling : " 
P1max=600 #MW
P1min=100 #MW
P2max=400 #MW
P2min=50 #MW
P3max=200 #MW
P3min=50 #MW
if P2<P2max and P2>P2min:
    print "P2 is within maximum and minimum limits." 
    P1=P1max #MW
    P3=P3min #MW
    P2=P-P1-P3 #MW

#Lambda=dC2/dP2 as P2 is niether maximum limit nor minimum limit.
dC2BydP2=7.8+2*0.0019*P2 #Rs./MWh
lamda=dC2BydP2 #Rs./MWh
dC1BydP1=6.5+2*0.0013*P1 #Rs./MWh
dC3BydP3=8.1+2*0.005*P3 #Rs./MWh
if dC1BydP1<lamda :
    print "Condition for P1 satisfied." 

if dC3BydP3>lamda:
    print "Condition for P3 satisfied." 

print "Load distribution is : " 
print "P1 = %0.2f MW" %P1
print "P2 = %0.2f MW" %P2
print "P3 = %0.2f MW" %P3

(b) According to optimum scheduling, Load distriution is :
P1 = 669.73 MW
P2 = 116.13 MW
P3 = 14.13 MW
(c) Optimum scheduling : 
P2 is within maximum and minimum limits.
Condition for P1 satisfied.
Condition for P3 satisfied.
Load distribution is : 
P1 = 600.00 MW
P2 = 150.00 MW
P3 = 50.00 MW

exa 2.18 page 62

In [23]:
import numpy as np
Bmn=np.mat([[0.0676, 0.00953, -0.00507],[0.00953 ,0.0521, 0.00901],[-0.00507, 0.00901, 0.0294]]) #Loss Coefficient
Bno=np.mat([[-0.0766],[0.00342],[0.0189]]) #Loss Coefficient
Boo=0.04357 #Loss Coefficient
P1=107.9 #MW
P2=50 #MW
P3=60 #MW
#solution : 
#PL=np.mat([[P1], [P2], [P3]])*Bmn+np.mat([[P1], [P2], [P3]])*Bno+Boo #MW
PL=np.mat([P1, P2, P3])*Bmn+np.mat([P1, P2, P3])*Bno+Boo #MW
pl= 0
for x in range(0,1):
    for y in range(0,3):
        pl+= PL[x,y]
print "Transmission Loss = %0.2f MW "%-pl 
#Note : Values calculated in the book are slightly wrong because of accuracy in calculation as compared to scilab accuracy.

Transmission Loss = 7.44 MW

exa 2.19 page 64

In [24]:
from sympy import symbols, solve
#lambda1=0.1*P1+20 #Rs./MWh
#lambda2=0.12*P2+16 #Rs./MWh
P=180 #MW
#Let loads are P1 & P-P1
#Economical loading lambda1=lambda2
P1, P2 = symbols('P1 P2')
eqn=0.1*P1+20-0.12*(P-P1)-16 
P1=solve(eqn, P1)[0] #MW
P2=P-P1 #MW
print "Load P1 = %0.2f MW "%P1 
print "Load P2 = %0.2f MW "%P2

Load P1 = 80.00 MW 
Load P2 = 100.00 MW

exa 2.20 page 65

In [25]:
from sympy import symbols, solve
#F1=0.004*P1**2+2*P1+80 #Rs./hr
#F2=0.006*P2**2+1.5*P2+100 #Rs./hr
P=250 #MW
P1, P2 = symbols('P1 P2')
dF1bydP1=2*0.004*P1+2 
dF2bydP2=2*0.006*P2+1.5 
#Let loads are P1 & P-P1
#Economical loading lambda1=lambda2
eqn=2*0.004*P1+2-2*0.006*(P-P1)-1.5 
P1=solve(eqn, P1)[0] #MW
P2=P-P1 #MW
print "Load P1 = %0.2f MW "%P1 
print "Load P2 = %0.2f MW "%P2

Load P1 = 125.00 MW 
Load P2 = 125.00 MW

exa 2.21 page 65

In [26]:
from sympy import symbols, solve
#F1=(8*P1+0.024*P1**2+80)*10**6 #Btu./hr
#F2=(6*P2+0.04*P2**2+120)*10**6 #Btu./hr
Pmax=100 #MW
Pmin=10 #MW
C=2.5 #Rs./million Btu
#C1=2.5*F1/10**6
#C2=2.5*F2/10**6
#For Maximum Load of 100 MW
P1, P2 = symbols('P1 P2')
dC1bydP1=8*2.5+2.5*2*0.024*P1 
dC2bydP2=6*2.5+2.5*2*0.04*P2 
#Let loads are P1 & Pmax-P1
#Economical loading lambda1=lambda2
eqn=8*2.5+2.5*2*0.024*P1-6*2.5-2.5*2*0.04*(Pmax-P1) 
P1=solve(eqn, P1)[0] #MW
P2=Pmax-P1 #MW
C1=2.5*((8*P1+0.024*P1**2+80)*10**6)/10**6 #Rs./hour
C2=2.5*((6*P2+0.04*P2**2+120)*10**6)/10**6 #Rs./hour
C100=(C1+C2)*12 #Rs.(Total cost of 12 hours on 100MW load)
#For Maximum load of 50 MW
#Let loads are P1 & Pmax-P1
#Economical loading : lambda1=lambda2
Pmax1=50 #MW
P1, P2 = symbols('P1 P2')
eqn=8*2.5+2.5*2*0.024*P1-6*2.5-2.5*2*0.04*(Pmax1-P1) 
P1=solve(eqn, P1)[0] #MW
P2=Pmax1-P1 #MW
C1=2.5*((8*P1+0.024*P1**2+80)*10**6)/10**6 #Rs./hour
C2=2.5*((6*P2+0.04*P2**2+120)*10**6)/10**6 #Rs./hour
C50=(C1+C2)*12 #Rs.(Total cost of 12 hours on 50MW load)
C=C100+C50 #Rs.(Total cost for 24 hours)
print "Minimum total cost for 24 hours = %0.2f Rs. "%C 
E=(Pmax*12+Pmax1*12)*10**3 #kWh
#Operating cost per unit energy
Co=C/E #Rs./kWh
print "Operating cost per unit energy = %0.2f Rs./kWh "%Co 
#Answer is wrong in the textbook. Calculation mistake in energy generation calculation & Cost calculation.

Minimum total cost for 24 hours = 49312.50 Rs. 
Operating cost per unit energy = 0.03 Rs./kWh

exa 2.22 page 66

In [27]:
from sympy import symbols, solve
import numpy as np
#F1=0.05*P1**2+21.5*P1+800 #Rs./hr
#F2=0.1*P2**2+27*P2+500 #Rs./hr
#F3=0.07*P3**2+16*P3+900 #Rs./hr
PT=200 #MW
Pmax=120 #MW
Pmin=39 #MW
#coefficients : 
c1=0.05; c2=0.1; c3=0.07 
b1=21.5 ;b2=27 ;b3=16 
a1=800 ;a2=500 ;a3=900 
lamda=(1/2*(b1/c1+b2/c2+b3/c3)+PT)/(1/2*(1/c1+1/c2+1/c3)) 
#Economical loading dF1/dP1=dF2/dP2=dF3/dP3
#from sympy import symbols, solve
P1=symbols('P1') 
P2=symbols('P2') 
P3=symbols('P3') 
dF1bydP1=2*0.05*P1+21.5 
dF2bydP2=2*0.1*P2+27 
dF2bydP3=2*0.07*P3+16 
#Solving equation :
A=np.mat([[2*0.05 ,0 ,0] ,[0 ,2*0.1, 0], [0, 0, 2*0.07]]) 
B=np.mat([[lamda-21.5], [lamda-27], [lamda-16]]) 
X=A**-1*B 
P1=X[0] #MW
P2=X[1] #MW
P3=X[2] #MW
if P2<Pmin:
    P2=Pmin 

P1plusP3=PT-P2 #MW
#dF1/dP1=dF3/dP3
#Let loads are P1 & P1plusP3-P1
P1=symbols('P1') 
P3=symbols('P3') 
eqn=2*0.05*P1+21.5-2*0.07*(P1plusP3-P1)-16 
P1=solve(eqn, P1)[0] #MW
P3=P1plusP3-P1 #MW
print "Optimum scheduling :" 
print "Loads P1, P2 & P3 are %0.2f, %0.2f & %0.2f MWs " %(P1, P2, P3)
F1=0.05*P1**2+21.5*P1+800 #Rs./hr
F2=0.1*P2**2+27*P2+500 #Rs./hr
F3=0.07*P3**2+16*P3+900 #Rs./hr
C=F1+F2+F3 #Rs/hour
print "For this schedule, total cost per hour = %0.2f Rs./hour " %C

Optimum scheduling :
Loads P1, P2 & P3 are 71.00, 39.00 & 90.00 MWs 
For this schedule, total cost per hour = 7190.65 Rs./hour

exa 2.23 page 67

In [28]:
from sympy import symbols, solve
#dF1/dP1=0.025*P1+15 #
#dF2/dP2=0.05*P2+20 #
PL=15.625 #MW
P1=125 #MW
lamda=24 #Rs. per MWh
B11=PL/P1**2 #Coefficient Loss
#dF2/dP2*L2=lambda
P2=symbols('P2') 
L2=1 #penalty factor
eqn=(0.05*P2+20)*L2-lamda 
P2=solve(eqn, P2)[0] #MW
#PL=B11*P1**2
P1=symbols('P1') 
dPLbydP1=2*B11*P1 
L1=1/(1-dPLbydP1) #penalty factor
eqn=(0.025*P1+15)-lamda/L1
P1=solve((eqn), P1)[0] #MW
print "Generation P1 & P2 are %0.2f & %0.2f MW"%(P1, P2)
PL=B11*P1**2 #MW
LD=P1-PL+P2 #MW
print "Load Demand = %0.2f MW "%LD 

Generation P1 & P2 are 123.29 & 80.00 MW
Load Demand = 188.09 MW

exa 2.24 page 68

In [29]:
from sympy import symbols, solve
#dC1/dP1=0.02*P1+16 #
#dC2/dP2=0.04*P2+20 #
PL=10 #MW
P1=100 #MW
lamda=25 #Rs. per MWh
B11=PL/P1**2; B22=0; B12=0 #Coefficient Loss
#dF2/dP2*L2=lambda
P2=symbols('P2') 
L2=1 #penalty factor
eqn=(0.04*P2+20)*L2-lamda 
P2=solve(eqn, P2)[0] #MW
#PL=B11*P1**2
P1=symbols('P1') 
dPLbydP1=2*B11*P1 
L1=1/(1-dPLbydP1) #penalty factor
eqn=(0.02*P1+16)-lamda/L1
P1=solve((eqn), P1)[0] #MW
print "Generation P1 & P2 are %0.2f & %0.2f MW"%(P1, P2)
PL=B11*P1**2 #MW
LD=P1-PL+P2 #MW
print "Load Demand = %0.2f MW "%LD

Generation P1 & P2 are 128.57 & 125.00 MW
Load Demand = 237.04 MW