Kothari D.P., Nagrath I.J. Modern Power Systems Analysis

G,5

CONTROL

STRUCTURES

IF

Statement

The

general

from

of

the

IF

statement

is

IF

expression

statements

ELSE

expression

statements

ELSEIF

statements

END

Expression

is

a logical

expression

resulting

in

an

answer,true,(l)

or

'false'(0).

The

logical

expression

can

consist

of

(i)

an

expression

containing

relational

operators

tabulated

along

with

their

meanings

in

Table

G.l.

Table

G.1

Relational

Operator

Meaning

Greater

than

Greater

than

or

Equal

to

Less

than

Less

than

or

Equal

to

Equal

to

Not

equal

to

(ii)

or

many

logical

expressions

combincd

with

Logical

operators.

Various

logical

operators

and

their

meanings

are given

in

Table

G.2.

Table

G.2

Logical

Operator

Meaning

AND

OR

NOT

FOR

Loops

This

repeats

a block

of

statements

precletermined

number

of times.

The

most

common

{&m

of

FOR

loop

used

is

fork=a:btc,

statements

end

where

k

is

the

loop

variable

which

is

initialised

to

value

of

initial

variable

a.

If

the

final

value (i'e.

c)

is

not

reached,

the

statements

in

the

body

for

the

loop

&

t -'__

I 041/

comes

to an end

when k reaches

or exceeds

the

final value

c.

For

example,

for

i

-

1:1:10,

a(i)

-

1

end

This

initializes

every

element

of a to

1. If increment

is

of 1,

as in

this

case,

then

the

increment

part

may as well

be

omitted

and

the above

loop

could

be

written

AS

for

i

-

1:10,

a(i)

=

1

end

While Loop

This

loop repeats

a block of

statements

till the

condition given

in the

loop

is

true

while

expression

statements

end

For example,

j

-1

while

i

<=

10

a(i)

=

1

i

=

i

+

1;

end

'

This

loop makes

first ten

elements

of array

a equal

to

l.

break statement

This

staternent allows

one to exit

prerr'aturely

from

a for

or while

loop.

G.6 HOW

TO

RUN THE PROGRAMS

GIVEN

IN

THIS

APPENDIX?

1. Copy these programs

into

the work

subfolder

under

MATLAB

foider.

2.

Just type

the name

of

the program

without'.m'

extension

and the program

will run.

3. If

you

wanL

to copy

them

in some

other

folder

say

c:\power,

then

after

copying those

files

in c:\power.

change

the work

folder to

e

:\power.

you

can

do this

by clicking on

toolbar

containing

three

periods

...

which

is

on

the right

side

to the Current

Directory

on the

top

rilht

corner.

4.

You can see

or edit

these

programs

by

going

through

Fite

-

open

menu

and

opening

the appropriate

file. However

do

not save

those programs,

unless you

are

sure that

you

want

the

changes

you

have

made

to

these

orisinal files.

ffiffi|

Modern

Power

system

Analysis

5. You can see which are the variables already defined

by

typing

whos

in

front of command

prompt.

That is why

you

will normally find a clear

command at the beginning

of our

programs.

This clears all the

variables

defined

so far from the memory, so that those variables do not interfere

G.7 SIMULINK BASICS

SIMULINK is a software

package

developed

by MathWorks

Inc.

which is one

of the most widely

used software

in academia

and industry for modeling and

simulating dynamical

systems. It can be used for modeling linear and nonlinear

systems, either

in continuous time frame or sampled time frame or even a hybrid

of

the two. It

provides

a

very

easy drag-drop type

Graphical

user

interface

to

build

the models in block diagram form. It has many

built-in block-library

components

that

you

can

use to model complex systems.

If these

built-in

models

are

not enough for

you,

SIMULINK allows

you

to have user defined blocks as

well. However, in this short appendix,

we will try to cover some of the

very

common

blocks that one comes across while simulating

a system.

You can try

to construct the

models

given

in the examples.

How

to Start?

You can start

SIMULINK

by

simply

clicking

the simulink icon in the tools bar

or by typing

Simulink

in front of the

MATLAB

command

prompt

>>.

This

opens up SIMULINK'hbrary browser, which should look similar

to

the

one

shown in Fig.

G2.

There may be other tool boxes

depending

upon

the license

you

have, The

plus

sign

that

you

see in

the

right half of

the

window indicates

that there

are

more blocks available

under

the icon clicking on the

(+)

sign will

expand the library.

Now for building up

a

new model

click on.File and select

New Model. A blank model window is opened. Now all

you

have to do is

to

select the

block

in the SIMULINK

library

browser and drop

it

on

your

model

window.

Then

connect them

together and run the simulation. That is

all.

An Example

Let us try to simulate a simple

model

where we take a sinusoidal input,

integrate

it and

observe

the

output.

The

steps

are

outlined as below.

1.

Click

on the Sources

in

the SIMULINK

library

browser window.

2. You are able to see

various

sources that SIMULINK

provides.

Scroll

down and

you

will see a Sine Wave sources icon.

3. Click on this sources

icon and without

releasing the mouse button

drag

and drop it in

your

model window which is currently named as

'untitled'.

4. If

you

double-click on

this

source,

you

will be able

to

see

Block

parameters

for sine wave which includes amplitude, frequency,

phase,

etc.

Let not change these

parameters

right now. So click on cancel to

go

back.

Fig.

G.2

Similarly click on

continuous library

icon.

You can

now

see various

built-in

blocks such as derivative,

integrator,

transfer-function,

state-

space etc.

Select

integrator

block

and drag-drop

it

in

your

model window.

Now click on sinks and drag-drop

scope

block

into

your

model.

This

is

one of the most common blocks

used for

displaying

the

values

of

the

blocks.

7.

Now

join

output of

sine-wave source to input

of

Integrator

bloc(.

This can

'

be done in two ways.

Either

you

click

the left

button and

drag

mouse from

output

of sine-wave source to input of

Integrator

block

and

leave left

button or otherwise click

on

right

button

and

drag the

mouse

to

form

connection

from input of integrator

block to

output

of sine,wave

source.

8. Now in the main menu, click

on

Simulation

and click

Start. The

simulation runs and stops after the

time specified

by

giving

ready

prompt

at

bottom left'corner.

9.

Now double-click on

scope

to

see the

output.

Is something

wrong?

The

result

is a sine wave of magnitude 2.

Is there

something

wrong with

SIMULINK software?

No, we have in fact forgotten

to

specify

the integration

constAnt!

Integration

of sin ?is

-

cos d+

C.

At

0=0,

C

-

1. If we

do

not specify

any

initial condition for output of

the integrator,

simulink

assumes

it to

be

0 and calculates the constant.

So it

calculates

flH

i

U

ccrlrnn

i".Sl

oircrctc

'.

S

furtmsaroUos

Unu

,

H

nmr'oer

i

.fi

sgn*aryrtan

iHs*'

:.

fl

sar:er

f,l

crrtrd 9y*en tooDox

Dbcret

Fqdin

t

l$hr

Ma['

Nmlle

Sign*

t SJ,s{emt

Si*r

1

:l

5.

6.

ffi

Modern Po

-

cos 0+

C

=

0

at

r

=

0

givinl

C

=

1. So the equation for output becomes

-

cos 0 +

l. Thus naturally, it starts from

0 at t

=

0 and

reaches its

peak

value

of 2 at

0= n,

i.e.

3.I4.

10. To

rectifv

this error. double click on Inte block and in the initial

conditions

enter

-

1 which should

be the output of the block at t

=

0. Now

run the simulation

again and see for

yourself

that

the

result is correct.

Some

Commonly

used

Blocks

1. Integrator

We have already described

its use in the above example.

2.

Transfer

function

Using

this block

you

can simulate a transfer function

of the

form Z(s)

=

N(s)/D(s), where

N(s) and D(s) are

polynomials

in s.

You

can double click the

block and enter the coefficients of s in numerator

and denominator

of the expression in ascending

order

of

s

which are to

be

enclosed

inside

square brackets and separated by.a

space.

3. Sum You can find this block under

Math

in

Simulink

block. By default,

it has two

inputs with both

plus

signs.

You

can

modify it to have required

number

of inputs to be

surnmed up by specifying a string of + or

-

depending

upon the inputs. So if there are

3

inputs

you

can

give

the list

of signs as +

- -

. This

will

denote

one

positive

input and two inputs

with

-

signs which

are often used to simuiate negative

feedback.

4.

Gain This block is also found under

Math in Simdlink block. It is used

to

simulate

static

eain.

It can even have fractional

values to act

as

attenuator.

5.

Switch This block is available under

Non-linear block in Simulink. It

has

3

inputs

with the top irrput being numbered

1. When the input number

2 equals or exceeds the threshold value

specified in the

properties

of this

block, it allows input number 1 to

pass

through, else it allows input

number

3 to

pass

through.

6. Mux

and Demux These blocks are available under

Signals and systems

block in

Simulink.

The Mux block combines

its inputs into a single

output

and

is mostly used to form a vector

out of

input

scalar

quantities.

Demux

block

does the reverse thing. It splits

the vector

quantity

into multiple

scalar outputs.

7. Scope We have already described its

use. However; if

you

are

plotting

a

large number of

points,

click on

properties

toolbar

and select Data

History

tab. Then uncheck the Limit

data

points

box so that all

points

are

plotted.

Also when the rvaveform does

not appear smooth, in the

general

tab

of

properties

toolbar, select

sample time instead

of

decimation

in the

sample

time box

and

enter a suitable value like

le-3

(10-3).

The scope

then uses the value at sample-time interval

to

plot.

8.

Clock This block is used to supply

time as a source input and is

available

in Sources block under

Simulink.

,*rr-rr r^t

in Sources block under

Simulink. We have

made

use of this

block in

stability studies to

provide

constant mechanical

input.

generate

a sine

wave

of any amplitude, frequency

and

phase.

The

reader is encouraged

to

work

out, the

examples given

in this

Appendix to

gain greater

insight into

the software.

G.8 SCRIPT AND FUNCTION FILES

Types of m-files

There are two types

of

m-files

used in Matlab

programming:

(i)

Script

m-file

-

This

file needs no input

parameters

and

does

not

return

any values as output

parameters.

It is

just

a set

of Matlab

statements

which

is

stored in a file so that one can

execute this

set by

just

typing the

file name in front of the command

prompt,

eg.

prograrnmes

in

G2

to Gl8

are script files.

(ii)

Function m-files

-

This file

accepts input agruments

and

return

values as

output

parameters.

It

can

work

with

variables

which

belong to the

workspace

as well as

with

the

variables which

are local

to the

functions.

These are useful for making

your

own function

for

a\

particular

application, eg. PolarTorect.m in

Gl

is a

function

m-file.

The basic

structure

of

function m-file is

given

below:

(a)

Function definition line

-

This is

the first line

of a function.

It

specifies

function name,

number and

order of input variables.

Its syntax

is- function

[output

variables]

=

function-name

(list

of input

variables)

(b)

First line of help

-

Whenever help is requested

on this funciton

or look

for is executed MATLAB displays

this

line.

Its

syntax is

-

Vo

function-name

help

(c)

Help

text

-

Whenever help is requested

on the

function-name

help

text

is displayed by Matlab in addition to first

help line.

Its syntax is

-

Vo

function-name

(input

variables)

(d)

Body

of the function

-

This consists

of codes

acting

on the set

of

input

variables to

produce

the output variables.

Tlccr nqn fhrrs rlcwelnn hic/hcr nrrrn rtrntrrqrnc end fi'rnnfinnc

qnA qAA

fham fn fha

LarvrII

uv ulv

existing library of functions and blocks.

G.9

SOME SAMPLE EXAMPLES

SOLVED BY

MATLAB

In

this

section 18 solved examples

of

this

book

are solved

again using

MATLAB/SIMULINK

to encourage the

reader

to solve more

power

system

problems

using MATLAB.

ffil

uodern

power

Systefn nnalysis

Ex

G.I

%

This

function

Yo

The

argument

% This

function

converts

po1

ar

to

this

function

has

been used i

to

rectangular

coordinates

is

1) Magnitude

& ?)Angle

is

.in

degrees

appendi

x

function

rect=polarTorect (a,b)

rect=a*cos

(b*pi/1S0)+j*a*si

n

(b*pi

/lg0);

function

is

used

in

solution

of many

of

the

examples

solved

later.

Ex

G.2

(Example

S.O)

% Thls

illustrates

the Ferranti

effect

% It

simulates

the

effect

by varying

the length

of

line from

%

zero(receiving

end)

to

5000km

in

steps

of

10

km

%

and

plots

the

sending

end voltage phasor

%

This

corresponds

to Fig.

5.13

and

data

from

Example

5.6

cl

ear

VR=220e3/sqrt

(3)

;

a1

pha=0.163e-3;

beta=1

.068e-3;

I

=5000;

k=1;

for

i=0:10:1,

VS=(VR/z)*exp(a'lpha*i)*exp(j*beta*j)+(VR/Z)*exp(-alpha*i)*exp(j*beta*j);

X(k)

=real

(VS)

;

Y(k)=imag(VS);

'

k=k+1;

end

p1

ot

(X,

Y)

Ex.

G.3

(Example

8.7)

% This

Program

illustrates

the

use

of

different

line

models

as

% in

Examp'le

5.7

cl

ear

f=50

|

=300

7=4Q+j*lZ5

Y=

I e-3

PR=50e6/3

YR=?20e3/

(sqrt

(3)

)

pf1

oad=0.8

IR=PR/(VR*pf1

oad)

z=7/1

y=Y

l1

%

Now we

calculate

the

sending-end

voltage,

send.ing-end

current

%

and

sending-end pf.by

following

methods

for

various

lengths

oF Iine

'

Appendix

G

ffi

T-

% varying

line lengths

from

10

to

300 km in

steps

of 10

km

% and

compare them

as a function

of line

lengths.

% The

methods

are

% l)

Short Line

approximatiom

Nominal-PI method

%

3) Exact transmission

line

equations

% 4)

Approximation

of

exact equations

i

=1.

for

I

=10:10:300,

%

Short

Iine

approximation

Vs_shortl'i ne

(i

)

=Yq+

(2*1

)

*I

R

Is shortline(i)=1P

spf_shortl i ne

(

i

)

=cos

(ang1

e

(vs_shortl

i ne

(i

)

-angl

e

(

I

s_shortl

ine

(

i

) ) ) )

Spower_shortl'i

ne

(i

)

=real

(Vs

shortl i ne

(

j

)

*conj

(

Is_shortl

i ne

(

i

)

) )

%Nomi nal PI

method

4=1+

(y*1

)*

(z*1)

/Z

D=A

B=z*l

C=y*t

*

(

l+

(y*1

)*

(z*1)

/

a)

Vs_nomi nal

pi

(

i

)

=A*VR+B*I

R

I

s_nom'i na'l

pi

(i

)

=C*VR+D*l

R

Spf_nominalpi(i)=cos(ang1e(Vs_nominalp'i(i)-angte(Is_nominalpi(i))))

Spower_nomi na'l

pi

(

i

)

=real

(Vs_nomi

nal

pi

(

i

)

*conj

(

I s_norpi

na I

pN.i

)

) )

% Exact transmjssion

Line Equatjons

Lc=sqrt(z/y)

garnma=sqrt (y*z)

Vs_exact(i)=cosh(gamma*l

)*VR

+

Zc*sinh

(ganrma*1

)*IR

Is_exact

(i)=(Illc)*sinh(gamma*1

)

+

cosh(gamma*l

)*IR

Spf_exact

(

i

)

=cos

(ang1

e

(Vs_exact

(

i

)

-ang1

e

(

I s_exact

(

j

)

) ) )

Spower_exact

(i

)

=real

(Vs_eiact

(i

)

*conj

(Is_exlct

(i

) ) )

% Approximation

of above

exact equations

4=1+(y*l

)*(z*t)

/Z

D=A

g=(z*1

)*(t+(y*l

)*(z*t

)/6)

C=

(y*l

)*(1+(y*l

)*

(z*1)

/6)

Vs_aPProx

(i

)

=A*VR+B*IR

I

s_approx

(

i

)

=[*!P+P*

1 P

Spf_approx

(i

)

=cos

(ang1

e

(Vs_approx

(

I

)

-ang

I e

(

I

s_approx

(i

) ) ) )

Spower_appnox(i

)

=rea1

(Vs_approx

('i

)*conj

(Is_approx(i

) ) )

point(j)=i

i=i+L

end

% The

reader

can uncomment

any

of

the

four

plot

statements

given

below

%

by removing the

percenlage

sign against

that

statement

% for

ex. in the

plot

statement

uncommented

below

% it

plots

the

sending end vo'ltages

for

short line

model

in red

%

by

nominal

pi-mode1

in

green

,

by

exact parameters

moder

%

and

by

approx.

pi

model

in

blue

pl

ot

(po'int,

abs

(Vs_shortl

i

ne)

,

'

r'

,

poi

nt,

abs

(Vs_nomi

na.lpi

)

,

(Vs_exact),

'b',

D

poi

nt,

abs

(Vs_approx)

,

k

,)

,

poi

nt,

abs...

,

poi

nt,

abs...

Ex

G.4

(Eiample

6,2)

%

This program

forms

yBUS

by

Singular

Transformation

&

The

data

for

this

program

is

a

primitive

admittance

matrix

y

%

which

is

to

be

given

in

the

foilowing

format

and

stored

in

ydata.

%

Ground

is

given

as

bus

no

0.

%

If

the

element

is

not

mutual]y

coupled

with

any

other

element,

%

the('the

entry

correspondi

ng

to

4th

and

5th

coi

umn

of

ydata

%

has'Io

be

zero

in

black

'

,

poi

nt,

abs...

%

element

no

I

connected

I

y (self)

lMutually

I

v(mutual)

?

|

From

I

To

I

tcoupled

tol

%

|

Busno

I

Busno

I

%

ydata=[

I

2

3

4

5

%

form

primitive

%

to

start

with

t

2

L/(0.05+;*9.15;

I

3

t/(0.1

+j*0.3

)

2

3

t/

(0.1s+j*0.4s)

?

4

1/

(0.10+j*9.39;

3

4

1/

(O

.05+j

*0.

15)

y

matrix

from

this

data

and

initialize

it

to

zero.

0

ol;

0

%

yo

%

el

ements=max (ydata

(

:

,

1)

)

ypri

mi

t i

ve=zeros

(e1

ements

,

e1

ements

)

%

Process

ydata

matr.ix

rowwjse

to

form

el

ement

no

I

connected

I

lFromlTo

I

Busno

I

eusno

I

y (sel

f)

1

2

L/

(O

.Os+j

*0.

15)

0

1

3

r/(o.i+j*0.3)

0

2

3

t/(0.r5+j*0.45)

0

2

4

t

/

(o

.19+3

*0

.30)

0

?, A 1 lln nc':+^ 1r\

v

r

r/\v.s3r..;"U.IC,f

U

lMutual

ly

I

y(mutual

)

I

coupl

ed

to

I

0

oj;

yprimi

ti ve

coupl

ed

wi

th

any

other

el

ement

i n

5th

col

umn

of

ydata

above)

the

i

th

row

i

s

made

equal

to

y(mutua.l

)

no

with

which

jth

element

is

mutually

coupled

for

i=1:elements,

ypri

mi

ti

ve

(i

,

.i)

=ydata

(i

,4)

%

Also

'if

the

element

is

mutually

eo

(whose

element

no

is

indicated

%

the

corresponding

column

no

in

t5

ydata=[

1

2

3

4

5

if

(ydata(i,S)

-=

g

;

j

is

ihe

element

j=ydata

(j

,5)

Ymutuai

=;'data(i,6)

yprimitive(i,j)

=

ymutual

end

% Form

Bus

'incidence

matrix

A

from

ydata

% Fol1owjng

statement gives

maximum

no.

of

elements

by

calculating

the

% maximum

out

of

all

elements

(denoted

by

:) of

first

column

of

yaata

el ements=max

(ydata

(:

,

l)

)

\

%

this

gives

no.

of

buses

which

is

nothing

but

the

maximum

entry

%

out

of

?nd

and

3rd

column

of

ydata

which

is

'from'

and

,to'

aolrrn,

buses=max

(max

(ydata

(z

,Z)),max

(ydata

(

:

,3)

) ) ;

ffi

__f--

buses=max(max(ydata

(2,

:)),max

(ydata

(3,

:

) )

)

A=zeros

(e'l

emen

t

s

,

bu

ses

)

for

i

=1:el

ements

%

ydata(i,?)

gives

the

,from'

bus

no.

%

fhe

entry

corrosponding

to

column

corrosponding

to

this

bus

%

in

A

matrix

is

made

1

if

this

is

not

ground

bui

if

ydata(i,21

-=

g

A(i,ydata(i,Z;;=1.

end

%

ydata(i,3)

gives

the

'to'

bus

no.

% The

entry

corrosponding

to

column

corresponding

to

this

bus

% in

A

matrix

is

made

1

if

this

is

not

ground

bus

i f

ydata

('i

,3;-=

9

A(i,ydata(i,3;;=-i

end

YBUS=A'*yprimitive*A

Ex

G.5

(data

same

as

Example

6.2)

This

program

forms

yBUS

by,adding

one

erement

at

a

t.ime

The

data

for

this

program

is

a

primitive

admittance

matrix y

which

is

to

be

given

in

the

following

format

and

stored

in

ydata

Ground

is

given

as

bus

no.

0

If

the

element

is

not

mutually

coupled

with

any

other

element,

then

the

entry

corresponding

to

4th

and

5th

coiumn

of

ydata

has

to

be

zero

The

data

must

be

arranged

in

ascending

order

of

element

no.

ffi,W

uooern Power

System nnavsis

YBUS=zeros

(buses,

buses)

;

for

row=1:elements,

% if

ydata(row,5)

is zero

that means the

corresponding element

% is not mutually

coupled with any other

element

nppencix o

Hffi

-

YBUS(k,l)

=

YBUS(k,l)

-

ysm(2,2);

YBUS(l,k)

=

YBUS(k,l);

YBUS

(i

1

,

k)

=

YBUS

(i

1,

k)

YBUS

(k,

i 1)

=

YBUS

(i

1, k)

;

YBUS(II,I)

=

YBtll(j!,L)

YBUS

(l

,j

1)

=

YBUS

(j

1, I

)

;

YBUS(i1,.|)

=

YBUS(i1,1)

YBUS

(.|

,

i 1)

=

YBUS

(j

1,.|

)

;

YBUS(j1,k)

=

YBUS(j1,k)

YBUS(k,jl)

=

YBUS(j1,k)

;

end

if

i1

==

0 &

j1

-=0

& k-=0 &

l-=0

YBUS(j1,j1)

=

YBUS(i1,i1)

+

ysm(1,1);

YBUS(k,k)

=

YBUS(k,k)

+

ysm(2,2);

YBUS(l

,l)

=

YBUS(1,.|)

+

ysn(?,2);

YBUS(k,l)

=

YBUS(k,l)

-

ysm(2,2);

YBUS(l

,k)

=

YBUS(k,l

)

;

YBUS(j1,.|)

=

YBUS(i1,1)

+

ysm(1,2);

YBUS(.l,j1)

=

YBUS(il,l)

;

YBUS(j1,k)

=

YBUS(i1,k)

-

ysm(1,2)

;

YBUS(k,j1)

=

YBUS(i1,k)

;

end

if i1-=0 &

jl==Q

&

k-=0 & l-=0

YBUS(i1,

j1)

=

YBUS(i1,'i1)

+

ysm(1,1);

YBUS(k,k)

=

YBUS(k,k)

+

ysn(2,2);

YBUS(l,l)

=

YBUS(1,'l)

+

ysn(z,2);

YBUS(k,l)

=

YBUS(k,l)

-

ysm(2,2);

YBUS(l,k)

=

YBUS(k,l);

YBUS(i1,k)

=

YBUS(i1,k)

+

ysm(1,2);

YBUS(k,i

1)

=

YBUS('il,k)

;

YBUS(i1,.|)

=

YBUS(i1,1)

-

ysm(1,2);

YBUS

(l

,'i

1)

=

YBUS

(i

1,.|

)

;

end

jf

i1-=0

&

j1-=Q

&

k==0

&

l-=0

YBUS

(i

1, i

1)

=

YBUS

(i

1, i 1)

+

ysm(1,l.)

;

YBUS(j1,i1)

=

YBUS(i1,i1)

+ ysm(1,1);

YBUS(l,l)

=

YBUS(1,.|)

+

ysm(2,2);

YBUS

(i

1,

j

1)

=

YBUS

(i

1,i 1)

-

ysm(1,l.)

;

YBUS(j1,j1)

=

YBUS(i1,i1)

;

YBUS(j1,1)

=

YBUS(j1,1)

+

ysm(I,2);

YBUS(.l,j1)

=

YBUS(i1,1)

;

YBUS(i1,.|)

=

YBUS(i1,.|)

-

ysm(1,2);

YBUS(l,jl)

=

YBUS(i1,1);

end

if i1-=0 &

jl-=Q

&

k-=0 & l==0

YBUS(i1,i1)

=

YBUS(i1,i1)

+

ysm(1,1);

YBUS(j1,jl)

=

YBUS(j1,i1)

+

ysm(1,1);

YBUS(k,k)

=

YBUS(k,k)

+

ysm(2,2);

YBUS(i1,j1)

=

YBUS(i1,j1)

-

ysm(1,1);

+

ysm(1,2);

+

vsm(1.2):

jf

ydata(row,5)

==

0

i 1

=ydata(row,2);

j1

=ydata(row,3);

if i1

-=

0 &

jl

-=

0

YBUS(i1,'i1)

=

YBUS(j1,j1)

+

ydata(row,4);

YBUS(i1,jl)

=

YBUS(i1,j1)

-

ydata(row,4);

YBUS(j1,i1)

=

YBUS(i1,jl);

YBUS(j1,j1)

=

YBUS(j1,j1)

+

ydata(row,4);

end

if il

==

0

&

jl

-=Q

YBUS(i1,j1)

=

YBUS(i1,i1)

+

ydata(row,4);

end

if

il

-=

0 &

jl

==Q

YBUS(j1,j1)

=

YBUS(j1,j1)

+

ydata(row,4);

end

eo

if

ydata(row,5)

is

NOT

zero

that means the

corresponding element

% is

mutually

coupled with

element

given

in

ydata(row,5)

if

ydata(row,5)

-=

g

i 1

=ydata(row,2);

j

1

=ydata(row,3)

;

%

mutualwith

gives

the

element

no

with which the

current element

% is mutually

coupled wjth k and 1

glve

the bus nos between

%

which

the mutually

coupled element is

connected

mutual *.i11=ydata

(

i 1,5)

;

;q=ydata

(mutual

wi th

,2) ;

1

=ydata

(mutual

wi th,3)

;

zsI=I/ydata(row,4);

zs?=L/ydata

(mutual

w'ith,4)

;

zm

=l/ydata(row,6);

756=

[zs

1 zm

zn zs?l;

ysm='i

nv

(zsm)

;

% Fol'lowing if

block

gives

the

necessary

modifications in YBUS

% when none of the buses is reference

(ground)

bus.

if il

-=

0 &

jl

-=

0 &

k

-=0

&

l-=0

YBUS('i1,i1)

=

YBUS(i

i,i

1)

+

ysm(1,1);

YBUS(j1,jl)

=

YBUS(j1,jl)

+

ysm(1,1);

YBUS(k,

k)

=

YBUS(k, k)

+

ysm(2,2);

YBUS(l,l)

=

YBUS(l

,l)

+

ysn(Z,2);

YBUS(i I,j1)

=

YBUS(i 1,j1)

-

ysm(1,1);

YBUS(j1,i1)

=

YBUS(i 1,j1)

;

-

ysm(1,2);

-

ysm(1,2);

W

Mooern

power

system

A_nalysis

YBUS(j1,i1)

= yBUS(it,jt);

YBUS(i1,k)

=

YBUS(i1,k)

+

ysm(1,2);

YBUS

(k,

i

1)

=

YBUS

(i

1,

k)

;

YBUS(i1,k)

= yBUS(i1,k)

_

ysm(1,2);

YBUS(k,il)

=

YBUS(j1,k)

;

end

YBUS

Ex

G.6

%

This'is

program

for

gauss

sieder

Load

frow.

The

data

is

cl

ear

n=4

V=[1.04

1.04

1

t]

Y=[3-i*9

-2+i*6

-l+5*3

0

-2+i*6

3666_j*11

_0.666+

j*2

_t+3*3

-1+;*3

-0.666+

j*2

3.666-

j*11

_Z+i*6

0

-1

+3'*3

-2+tr*6

3-j"91

type=ones

(n,1)

typechanged=zeros

(n,

1

)

Ql

i

mi

tmax=zeros

(n,

1

i

Ql

imi

tmi

n=zeros

(n,

l)

Vmagfi

xed=zeros

(n

,

1

)

tYpe(2)=2

Ql

i

mi

tmax

(2)

=1

.

g

Q1

i

mi

tmi

n

(2)

=9.

2

Vmagfi

xed(2)

=1.04

di

ff=10;

noofi

ter=1

Vprev=V;

while (diff>O.00001

|

noofiter==1),

abs

(V)

abs

(Vprev)

%pause

Vprev=V;

P=[inf

0.5

-1

0.3];

Q=[inf

0

0.5

-0.1]

;

$=[inf+j*inf

(0.5-j*0.2)

(-1.0

+

j*0.5)

(0.1-1*g.1)];

for

i

=Zin,

if

type(i)==2

ltypechanged(i;==1,

",llilil;:il;l#ii#ll,

ill

).Q,

imi

tmi

n(i

) ),

Q(j

)

=Ql

imi

tmi

n

(i

) ;

el

se

Q('i)=Qlimitmax(i);

end

tYPe

(i

)

=1;

from

Example

6.5

lFitlilE4lq.

ffi.

Appendix

c

NilFffi

T-

typechanged

(i

)

=1

t

el se

type

(i

)

=2;

end

sumyv=0;

for

k=L:n,

if(i

-=

k)

sumyv=sumyv+Y

(i,

k)

*V

(k)

;

end

v

(i

)

=

(1/Y

(i,

i

) )

*

(

(p (i

)

-j*Q(i

) )

/conj

(v

(i

)

-sumyv)

;

if

type(i

;==2

& typechanged(i)-=1,

V(i

)=pe1

arTorect(Vmagfixed(i

),ang1

e(V(i

)

)*1g0/pi

) ;

end

di ff=max

(abs

(U

(Zzn)

)

-abs

(Vprev

(Z:

n)

) ) ) ;

noofi ter=noofi

ter+l;

end

V

Ex.

G.7

(Example

6.6)

% Program

cl ear;

% n

stands

for

load

flow

by

Newton-Raphson

Method.

for

number

of

buses

n=3;

%

Y

,voltages

at

those

buses

are initjaljsed

v=[1.04

1.0

1.04];

% Y is

YBus

Y=[

5.88228-j't23.50514

-2.9427+j*Lt.t6t6

-2

.9427

+j*Lt

.7

67 6

5 .88228_j*23

.505L4

-?

.9427

+

j*Ll

.7

67 6

-2.9427

+

j*tt

.7 6t6

%Bus

types

are

jnitialjsed

in

type

array

to

%bus.

%code

2

stands

for

PV

bus

type=ones(n,1);

%when

Q

l'imjts

are

exceeded

for

a

pV

bus

Bus

type

is

changed

to

pQ

%temporari

1

y

Lan a'l omonf i nf #.r^oahrnan'l

.i

a -^+ .^ 1 i -

'e",L,,e

I

vr

uJpsLrsrvsu

r)

JtrL

Lu

r rlr

cdsg

rr5

Dus

status

%is temporarily

changed

from

pQ

to

pV.Otherwise

it

is

zero

typechanged=zeros

(n,

1)

;

%

since max

and

min

Q

l'imits

are

checked

only

for

pv

buses,

% max

& min

Q

limlts

for

other

types

of

buses

can

be

set

to

any

values.

% here

we

have

set

them to

zeros

for

convenience

Ql

imi

tmax=zeros

(n,

1)

;

Q1

imi

tmi n=zeros

(n,

l);

Vmagfi

xed=zeros

(n,

1)

;

-Z

.9427

+

j*

7I

.t

67

6

-2

.9427

+

j*

II

.7

67

6

5.88228-j*23.5051a1;

code

1 which

stands

for

PQ

W

Modern

Power system Analysis

I

t I to

total

no.

of equations

(eqcount)

for

ceq=1:eqcount,

for

ccol=1:eqcount,

vsJ\eey/tsJrvvv

r\uJJvevruuJ\vwwrllr,

bm=i mag

(Y (assoeqbus

(ceq),

assocol

bus

(cco:l

) )

*V

(as

socol

bus

(ccol

) ) )

;

ei

=real

(V (assoeqbus

(ceq)

) )

;

fi

=imag

(V

(assoeqbus (ceq)

) )

;

i f assoeqvar(ceq)

=='P

"'

& assocol vai(cco1

)

=='d'

,

i f assoeqbus

(ceq)

-=assocol

bus

(ccol

)

,

H=am*fi

-bm*ei

;

el se

11

=-Q

(as

soeqbus

(

ceq)

)

i mag

(Y (

assoeqbus

(ceq),

assocol bus

(ceq))

. . .

*abs

(V

(assoeqbus (ceq)

))^2) ;

end

Jacob

(ceq,

ccol

)

=H

end

i f assoeqvar(ceq)

=='P'

& assocol var(cco1

)

=='V'

,

'i

f assoeqbus

(ceq)

-=assocol

bus

(ccol

)

,

N=am*ei+bm*fi;

el se

N=P

(assoeqbus

(ceq)

)

+real

(Y (assoeqbus (ceq)

,

assocol bus

(ceq)

)

.

*abs

(V (assoeqbus (ceq)

) )^2) ;

end

Jacob

(ceq,

ccol

)

=11

end

if assoeqvar(ceq)=='Q'

& assocolvar(cco1)=='d',

i

f

assoeqbus

(ceq)

-=assocol

bus

(ccol

)

,

J=am*ei+bm*fi;

el se

J=P

(as

soeqbus

(

ceq

))real

(

Y

(assoeqbus

(ceq

),assocol

bus

(ceq)

)

*abs

(V (assoeqbus

(ceq)

) )^2)

;

end

Jacob(ceq,ccol

)

=J

end

if assoeqvar(ceq)=='Q'

&

assocolvar(cco1)=='V',

1 f assoeqbus

(ceq)

--assocol

bus

(ccol

)

,

L=am*fi

-bm*ei

'

el se

L=Q(assoeqbus

(ceq)

)

-.

. .

imrn/V/'cc^o^h,,.t/^o^\ 'e.^^^'1 h,,./^^^\\*.h./\, 1...^^^h,,-/^^^\\\'\D\.

rfffuv\r\.uJJvsqvuJ\lsY/rqJJvuvlpuJ\LsY/r/,

qwJ\r\qJJvEVUUJ\Lgr{,,

ll

Ll,

end

Jacob

(ceq,

ccol

)

=L

en0

end

en0

%New

Update vector is

calculated

from

Inverse of the

Jacobian

pause

update=i

nv

(Jacob)

*mi

smatch'

;

noofeq=L;

if

type(i;==1

newchi nangV=update

(noofeq)

;

newangV=ang1

e

(V (

i

) )

+newchi

nangV

;

newchi

nmagV=update(noofeQ+l)*abs

(V (i

) ) ;

newmagV=abs

(V

('i

)

+newch

j

nmagV;

V

(j

)

=pol

arTorect

(newmagV,

newangV*180/pi)

;

noofeq=noofeq+l

;

el

se

newchi

nangV=update

(noofeq)

;

newangV=ang1

e

(V ('i

) )

+newctri

nangV

;

V(i

)=pol

arTorect

(abs

(V(i

)),newangV*l8O/pi

)

;

noofeq=noofeq+l

;

end

% Al I

the

fol

I

owi

ng vari

abl

es/arrays

are

cl eared

from

% memory.

This

is because

the'ir dimensions

may

change

due

to

% bus

switched

and once updates

are

calculated,

the variables

9r

are of no use

as

they are being

reformulated

at the

% end of each iterat'ion

\

clear mismatch

Jacob update

assoeqvar

assoeqbus

assocolvaf

assocolbus;

di

f f=m'i n

(abs

(abs (V (2:

n)

)

-abs

(Vprev (2:

n)

) ) ) ;

noof

i ter=noofi

ter+1

;

end

Ex.

G.8

(Table

7.1)

% MATLAB

Program

for

opt'imum loading

of

generators

% The data

are from

Example

7.1

% rt finds

lamda by

the algorithm

given

on the

same

pd9€,

once

the

demand

rk

is

spec'ified

% We have taken

the

demand as

231.25MW

corrosponding

to

the

last

but

one

% row of Table

7.1 and calculated

lamda

and the

load

sharlng

%

n

is no

of

generators

n=2

% Pd

stands for

load demand.

%

alpha and beta

arrays

denote

alpha

beta

coeffjcients

%

for

given

generators.

Pd=?31.25

al

pha=

[0.20

0.251

beta=

[40

301

% initial

guess

for lamda

I

amda=20

m't n't mum

%are stored

in arrays

% In real

life large

% array

to inf

using

ffil

Modern

po

I

amdaprev=l amda

%

tolerance is

eps and increment

in lamda is

deltalamda

eps=1

de'l

tal amda=0.25

Hffi

m1n

"

lnfi

cost=0;

for

j=0:n,

for

i=0:

n,

for I

=0:n,

unit

=

[0

0 0 0];

%

Here

we elimnate

straightaway

those

combinations

which

% dont

make up

the

% n

MW demand

or

such

combinations

where maximum

generation

% on

lndividual

%

generation

is exceeding

the

maximum capacity of

any of the

%

generators

if('i+j+k+'l)==n

&

icPgmax(1)

&

i<Pgmax(2)

&

kcPgmax(3)

&

I

<

...

Pgmax

(4)

if

i-=0

unit(1,1)=i;

%

Ftnd out

the cost of

generating

these units and

%

qdd

'it

up to

total cost

cost=cost+0.

5*al

pha

(

I

)

*i *i

+beta (

I

)

*i

;

end

i f

3-=g

uni t

(

L,2)=i;

cost=cost+0.5*alpha(2)*j*j+beta(2)*i

;

\

end

i f k-=0

uni

t

(1,3)

=k;

cost=cost+0.5*al

phr(3)*k*k+beta(3)*k

;

end

i f l-=0

unit(1,4)=l;

cost=cost+0.5*a1

pha(a)*l

*l

+beta(4)*l

;

end

%

If

the total

cost

i

s comi

ng out to

be

I

ess

than

% minimum

of

the cost

in

%

then

make min equal

to cost

and

% cunit

(stand

for committed

units)

equal

to un'its

% committed

in

this

iteration

t

(denoted

bY

variable

units)

if cost

<

min

cun'it

-

un'iti

mi n=nnc.f :

rrr

. | | ev v

v

t

el se

cos

t

=0;

end

mum

|

'r

mt ts

ot each

generati

ng

un

Pgmin

and

Pgmax.

scale

problems,

we

can

first

initialse

the Pgmax

for

I

oop and

Lamda i s')

I

oad shared

by two un

j

ts

'is')

%

Pgnin

array

to zero using

Pgmin=zeros(n,1)

command

% Later we

can

change the

l'imits

indiv'idually

Pgmax=[125

I25]

Pgmin=120

Z0)

Pg--1O0*ones

(n,1)

whi

le abs

(sum(PS)

-ed;'gtt

for

i=1:n,

Ps(i

)

=

(1

amda-beta

(i

)

)/al

pha

(i

)

;

if

Pg(i

)>Psmax(i

)

Pg(i

)=Psmax(i

)

;

end

if Ps(i).Pgmln(i)

Pg(i)=Pgmin(i);

end

if

(sum(Pg)-Pd

).0

lamdaprev=l amda;

lamda=l amda+del ta1 amda;

el se

I amdaprev= I amda

;

lamda=l

amda-del tal amda;

end

disp(' The final value of

1

amdaprev

disp('

The

distribution of

Pg

Ex. G.9

(Table

7.2)

% MATLAB

Program

for optimum

unit

committment

by Brute Force method

ft The data

for

thls

program

corresponds

to

Table

7.2

c

I

ear;

% alpha

and beta

arrays

denote

alpha

beta coefficients for

given generators

alpha=10.77 1.60

2.00 2.501';

beta=123.5 26.5

30.0

32.01';

Pgmin=[1 1 I 1]';

Pgmax=ll? 12 12 I2l':'

n=9

% n denotes total MW to

be

commi tted

end

,f,f6.'l

Modern Po@

end

E

for i=1:n,

s

i

gma=B

(i

,

:

)

*Pg-B

(

i

,

i

)

*Pg

(

i

)

;

end

, disp

('cunit

display

the no

of

commjtted units

on each

of

tors')

disp('

If cunit

for a

particu)ar

generator

is

0

jt

means

for

I amda

the four

genera-

the unit is

not

comm

di

sp

(

'

The total

no

of units to

be commi

tted are'

)

cunit

Ex

G.lO

(Ex.

7.A)

clear

%

MATLAB

Program

for optimum

loading

of

generators

%

This

program

finds

the optimal loading

of

generators

including

%

penal

ty factors

% It

implements

the

algorithm

given

just

before Example

7.4.

%

The

data for this

program

are taken

from

Example

7.4

%

Here

we

give

demand

Pd and alpha,

beta

and B-coefficients

t

l,{e calculate

load

shared by

each

generator

%

n is

no of

generators

n=2

? Pd

stands for I

oad demand

%

alpha

and beta

arrays

denote alpha

beta

coefficjents for

given

%

generators

Pd=237

.04;

al

ph6=

[0.020

0.0a1;

beta=

[16

2ol

i'

I

initial

guess

I amda=20;

I amdaprev=

I

amda;

% tolerance

is eps

and increment

in lamda

js

deltalamda

eps=1;

del tal

amda=O.25;

%

the

minimum

and maximum

limits

of each

generating

unit

%are stored in

arrays Pgmin

and Pgmax.

% In actual large

scale

problems,

we

can first initialise

the Pgmax array

% to inf

using for

loop

% and Pgmin

array to

zero using

Pgmin=zeros(n,1)

command

% Later

we

should can

change the limits

indiv'idua11y

Pqmax=f200

200.l:

Pgmin=[0

0];

B=[0.001

0

0 0l;

Pg=zeros(n,1);

noof

i ter=0;

PL=0;

Pg=zeros(n,1);

Pg

(i

)

=

(1-

(beta (i

)/l

amda)

-

(2*s

i

sna))

/

(al

pha

(i

)

/l

amda+z*B

(i,

i

)

;

I

if Pg(i)'Pgmax(j)

Pg

(i

)

=Pgmax

(i

)

;

end

if Pg(i)<Pgmin(i)

Ps(i)=Psmin(i);

end

PL=Pg'*B*Pg;

if

(sum(Pg)-Pd-PL

)<0

1

amdaprev=l amda;

I amda=l amda+del tal amda;

el

se

I

amdaprev=1

amda;

I amda=l amda-del tal amda;

end

noofi ter=noofi ter+1

;

Pg;

end

di sp

(

'The

noofi ter

di

sp

(

'The

1 amdaprev

di

sp

(

'The

Pg

di

sp

(

'The

PL

Ex.

G.II

% Iri this

%

Fis.

8.8

s'imul ati on

Tsg=

.4

Tt=0.5

Tps=20

Kps=100

R=3

Ksg= 1g

Kt=0. I

Ki

=0.09

no

of iterations

required

are')

f i nal val ue

of

'lamda

i s

'

)

opt'ima1 I oadi

ng

of

generators

i ncl

udi ng

penal

ty facto,rs

'is'

)

I osses are')

example

the

parameters

are initialized This

both for Figs.

G.3 and

for

all the blocks

for

the system

in

program

has to

be run

prior

to

the

G.4.

Kothari D.P., Nagrath I.J. Modern Power Systems Analysis

Подождите немного. Документ загружается.