关于matpower中目标函数修改的问题

clear;

clc;

errArr=[];

%%

%³õʼ»¯£¡£¡£¡

initial;

% Start clock

t1 = clock;

%%

ROU=sl'*MU_MIN+su'*MU_MAX;

MUt=SIGMA*ROU/(2*length(sl));%³õʼ¶ÔżÒò×ÓÓë³Í·£Òò×Ó¼ÆËã%

ik=0;%¼Æµü´ú´ÎÊý£¡£¡£¡

%µü´úÑ »·¹ý³Ì£¡£¡

while(abs(ROU)>=err)

: S3 E& ~! p$ T' |' m. G
%%

0 c0 z! j5 `8 A" t: V3 {- x; @
%Calcute h,g matrix

    ROU=sl'*MU_MIN+su'*MU_MAX;

        errArr=[errArr;ROU;];

    SIGMA=0;

    MU=SIGMA*ROU/(2*length(sl));         %ÖÐÐIJÎÊýÖÃÁã%

$ e0 x/ K4 j0 ^( d
for i=1:30

        temp=0;

( f" U5 R9 o3 a( R
for j=1:30

            temp=temp-V(j)*aY(i,j)*cos(Vth(i)-Vth(j)-Yth(i,j));

6 N. |8 Y3 j$ j1 n3 yend

+ b- u* F" W+ F7 r. n# Cif (i>6)

            tPg=0;

; q4 ?7 S3 U) Yelse

            tPg=Pg(i);

: f) \% B' [  t6 j% ^
end

        h(i)=tPg-Pd(i)+V(i)*temp;

3 G4 f. M7 _: d# X/ U7 [8 _) `end

* S, x8 C+ t2 z, o
for i=1:30

        temp=0;

4 c5 [) H3 g% o9 V3 j- l
for j=1:30

            temp=temp-V(j)*aY(i,j)*sin(Vth(i)-Vth(j)-Yth(i,j));

& K# f; G+ q$ p+ C7 Z8 h& G5 a7 w5 q5 D
end

/ X$ O- x- M3 e8 I
if (i>6)

            tQg=0;

& }0 f- ], E- p/ }+ q# T/ r
else

            tQg=Qg(i);

' Q$ E+ Q/ m1 X3 `4 Q0 H% t: X. D
end

        h(i+30)=tQg-Qd(i)+V(i)*temp;

, M* _8 k# M( s8 @end
. k. P8 e( I* H0 f+ X! n% Cal  h END

; M, I) I/ Z0 W* U

4 l; _) K) |1 D9 ]& b, y
for i=1:6

        g(i)=Pg(i);

        g(i+6)=Qg(i);

' i* W. Z! J& s. w' ?; t. \4 {
end

) Q. t  L& {  K4 {* t" s  Y* b
for i=1:30

        g(i+12)=V(i);

8 f( T5 y- N( U' N, nend! e# h; |/ l7 A3 ^0 q: V
% Cal  g  END

! _6 V6 r; i  D$ A
%Calcute h,g matrix END

) F/ `% d' b# K/ y& W  ~
%%

) ]  t6 B! Z1 P) u, }, Y
%Calculate Jacobian&Hessian matix

1 u# g) {7 o. j8 i% C! H0 X1 ~' G$ h. o
%First Step: Jf,Hf

6 I  }8 J- v  j/ v6 }/ n5 B  a* Afor i=1:6

        Jf(i)=2*gencost(i,5)*Pg(i)+gencost(i,6);

        Hf(i,i)=2*gencost(i,6);

# f3 C6 K& O) F. L. T! p$ jend

, Y- Z" s1 l$ c5 a9 j2 @
%Second Step: Jh, hΪµÈʽԼÊø

  n7 ~- O& b5 v, k& ^! Q, K
for i=1:6         %Ç°6ÐжÔPgÇóµ¼£¬ÓÉ´ËÒÑÇó³ö

        Jh(i,i)=1;

, C+ h" K- [% B5 Z/ W5 P% R
end

" G' u/ x4 b9 @* D  Dfor i=7:12        %7-12ÐжÔQgÇóµ¼£¬ÓÉ´ËÒÑÇó³ö

        Jh(i,i+24)=1;

/ `6 ^; E6 b5 p' t
end

( h: T$ x6 J4 H  ~
for i=1:30       %ÐγÉ13-42ÐеÄ1-60ÁÐ

# D5 |5 w* C, _9 R$ r; Pfor j=1:30

            tempVp=0;

            tempVq=0;

2 V9 B" p/ v* U1 s' Y: S
if (j==i)

6 J, `/ {9 k% d- V9 K) k) K: Zfor k=1:30

                    tempVp=tempVp-V(k)*aY(j,k)*cos(Vth(j)-Vth(k)-Yth(j,k));

                    tempVq=tempVq-V(k)*aY(j,k)*sin(Vth(j)-Vth(k)-Yth(j,k));

1 i+ z+ B( z1 N9 z4 q
end

                Jh(12+j,i)=tempVp-aY(j,j)*V(j)*cos(Yth(j,j));

                Jh(12+j,30+i)=tempVq+aY(j,j)*V(j)*sin(Yth(j,j));

% C1 b, B  T/ p9 @1 w% E1 _2 Felse

                Jh(12+j,i)=-aY(i,j)*V(i)*cos(Vth(i)-Vth(j)-Yth(i,j));

                Jh(12+j,30+i)=-aY(i,j)*V(i)*sin(Vth(i)-Vth(j)-Yth(i,j));

; e4 s8 E3 G& t" _4 }end

! ^3 x7 e! e" V& v
end

6 g; s# Z0 J4 G3 H  X
end

# X) K; b  G8 T
for i=1:30       %ÐγÉ43-72ÐеÄ1-60ÁÐ

2 D1 S" {' W' i1 ~/ ?# f
for j=1:30

            tempVp=0;

            tempVq=0;

9 A! G8 h5 L; f/ y3 Cif (j==i)

+ I9 ~2 y3 V0 v8 W2 G& H
for k=1:30

                    tempVp=tempVp+aY(j,k)*V(k)*sin(Vth(j)-Vth(k)-Yth(j,k));

                    tempVq=tempVq-aY(j,k)*V(k)*cos(Vth(j)-Vth(k)-Yth(j,k));

( G: h9 h% B8 e  zend

                tempVp=tempVp-V(j)*aY(j,j)*sin(-Yth(j,j));

                tempVq=tempVq+V(j)*aY(j,j)*cos(-Yth(j,j));

                Jh(42+j,i)=V(i)*tempVp;

                Jh(42+j,30+i)=V(i)*tempVq;

  l4 u& {$ v2 q& G$ x
else

                Jh(42+j,i)=-aY(i,j)*V(i)*V(j)*sin(Vth(i)-Vth(j)-Yth(i,j));

                Jh(42+j,30+i)=aY(i,j)*V(i)*V(j)*cos(Vth(i)-Vth(j)-Yth(i,j));

4 P9 T/ [+ A" e6 L9 m) Oend

' k% `( n; Y/ K4 {2 m
end

/ q2 \/ z5 J. S+ v4 y
end

1 F  a% F% t; @# v2 q0 u- R
%Third Step: Hh

' g% V+ Z* I5 \6 Z. i
%Óй¦²¿·Ö

! C( t& ]$ M  `' p! e
for i=1:30

7 p8 }+ b2 L6 i2 y2 d5 V0 l* ifor j=1:30

$ m1 i6 r' f* F2 {$ C* @
for k=j:30

* ?: O$ U- Q, C# H4 _3 R
if (j==k&&i~=j)

                    Hh(j+12,k+12,i)=0;          %VV

                    Hh(j+42,k+42,i)=V(i)*aY(i,j)*V(j)*cos(Vth(i)-Vth(j)-Yth(i,j)); %%thth

1 q3 B9 |- B& a+ a- A! F
elseif (j==k&&i==j)

                    Hh(j+12,k+12,i)=-2*aY(j,j)*cos(Yth(i,i));  %VV

                    temp=0;                   %thth

, q! ~; l" p* o5 |for l=1:30

                        temp=temp+aY(j,l)*V(l)*cos(Vth(j)-Vth(l)-Yth(j,l));

2 ^, _: o" L5 t
end

                    temp=temp-aY(i,i)*V(i)*cos(-Yth(i,i));

                    Hh(j+42,k+42,i)=V(i)*temp;

+ A& c) X" }% p, M' Welseif (k==i)

                    Hh(j+12,k+12,i)=-aY(i,j)*cos(Vth(i)-Vth(j)-Yth(i,j));    %VV

                    Hh(k+12,j+12,i)=Hh(j+12,k+12,i);

                    Hh(j+42,k+42,i)=-V(i)*aY(i,j)*V(j)*cos(Vth(i)-Vth(j)-Yth(i,j));  %thth

                    Hh(k+42,j+42,i)=Hh(j+42,k+42,i);

3 k/ q; Y5 P$ Welseif (j==i)

                    Hh(j+12,k+12,i)=-aY(i,k)*cos(Vth(i)-Vth(k)-Yth(i,k));    %VV

                    Hh(k+12,j+12,i)=Hh(j+12,k+12,i);

                    Hh(j+42,k+42,i)=-V(i)*aY(i,k)*V(k)*cos(Vth(i)-Vth(k)-Yth(i,k));  %thth

                    Hh(k+42,j+42,i)=Hh(j+42,k+42,i);

& ~- T5 J) E5 |% T2 K$ xend

3 U; V" \( u# [0 m
end

* Y! a) ]# E' u; c) r
end

" R% n7 T( s9 v3 o' \5 n
end1 O( f+ r& P8 S' [5 R9 A
%ÖÁ´ËÒÑÐγɣ¨13-42£¬13-42£©ºÍ£¨42-72£¬43-72£©

2 B; Y1 c  ^9 e1 |6 `
for i=1:30

- l: t7 w9 N5 M7 ~5 ^+ \+ [# e$ Kfor j=1:30

3 k' _) o6 Z) K, M# n
for k=1:30

% j! C$ H5 M, G9 a
if (j==k&&i~=j)

                    Hh(j+42,k+12,i)=-V(i)*aY(i,j)*sin(Vth(i)-Vth(j)-Yth(i,j));      %thV

+ s2 U7 }0 v# C) q6 n8 y' relseif (j==k&&i==j)

                    temp=0;                   %thV

# V7 e% }% X$ X+ p6 ofor l=1:30

                        temp=temp+aY(j,l)*V(l)*sin(Vth(j)-Vth(l)-Yth(j,l));   

. a7 R7 n6 W! W0 ?, h
end

                    Hh(j+42,k+12,i)=temp-V(i)*aY(i,i)*sin(-Yth(i,i));

- l0 W4 b' R' j; x1 aelseif (j==i)

                    Hh(j+42,k+12,i)=V(i)*aY(i,k)*sin(Vth(i)-Vth(k)-Yth(i,k));   %thV   

- e% r' w3 @( ]& w
elseif (k==i)

                    Hh(j+42,k+12,i)=-V(j)*aY(i,j)*sin(Vth(i)-Vth(j)-Yth(i,j));  %thV

! c6 V/ e& ]9 o7 v/ U7 aend

: S) |9 i$ b$ N) j& zend

4 X* d' F7 q+ J+ ^2 U
end

        Hh(13:42,43:72,i)=Hh(43:72,13:42,i)';

" }8 V, q3 U# t0 ?1 |+ l3 z% N" R4 R
end
: D& O% D" M/ c%ÖÁ´ËÒÑÐγɣ¨42-72£¬13-42£©ºÍ£¨13-42£¬43-72£©

$ i: u4 S0 g! L; P  g
%ÎÞ¹¦²¿·Ö

2 c- Z0 r3 ]+ c: Q% r; \  `$ Z
for i=1:30

! G4 R0 x4 t, h9 Z! g) ]( jfor j=1:30

& F8 _) Z7 @2 Y+ |5 z0 ~; w0 O+ ~
for k=j:30

2 _% _& D: `" o% ~" M5 T3 F- H& }if (j==k&&i~=j)

                    Hh(j+12,k+12,i+30)=0;          %VV

                    Hh(j+42,k+42,i+30)=V(i)*aY(i,j)*V(j)*sin(Vth(i)-Vth(j)-Yth(i,j)); %%thth

$ z8 `) L; ~3 I* ^" v- ~elseif (j==k&&i==j)

                    Hh(j+12,k+12,i+30)=2*aY(j,j)*sin(Yth(i,i));  %VV

                    temp=0;                   %thth

0 j, y$ E/ L7 Y2 c$ C! l/ q" ofor l=1:30

                        temp=temp+aY(j,l)*V(l)*sin(Vth(j)-Vth(l)-Yth(j,l));

8 O+ n1 M, D* x, Qend

                    temp=temp-aY(i,i)*V(i)*sin(-Yth(i,i));

                    Hh(j+42,k+42,i+30)=V(i)*temp;

* w- `) j( S. N8 V5 j' Aelseif (k==i)

                    Hh(j+12,k+12,i+30)=-aY(i,j)*sin(Vth(i)-Vth(j)-Yth(i,j));    %VV

                    Hh(k+12,j+12,i+30)=Hh(j+12,k+12,i+30);

                    Hh(j+42,k+42,i)=-V(i)*aY(i,j)*V(j)*sin(Vth(i)-Vth(j)-Yth(i,j));  %thth

                    Hh(k+42,j+42,i+30)=Hh(j+42,k+42,i+30);

& o& @: `" _2 D6 h2 T, `* ]elseif (j==i)

                    Hh(j+12,k+12,i+30)=-aY(i,k)*sin(Vth(i)-Vth(k)-Yth(i,k));    %VV

                    Hh(k+12,j+12,i+30)=Hh(j+12,k+12,i+30);

                    Hh(j+42,k+42,i+30)=-V(i)*aY(i,k)*V(k)*sin(Vth(i)-Vth(k)-Yth(i,k));  %thth

                    Hh(k+42,j+42,i+30)=Hh(j+42,k+42,i+30);

: U/ d1 _) d( n4 gend

' u7 i3 s: r) i( h! G3 t+ p1 S
end

- f# S2 b+ P% ~1 Q3 P2 mend

* `0 G/ i$ d1 x. w: }end& D) E# I5 d" c5 a3 c! [/ i
%ÖÁ´ËÒÑÐγɣ¨13-42£¬13-42£©ºÍ£¨42-72£¬43-72£©

3 @9 g5 P8 v* i8 r. I' J/ U
for i=1:30

& P8 G( O; E% b9 z
for j=1:30

' s& @. Q. X' l6 S3 u3 Vfor k=1:30

" J7 z7 d7 b+ l5 `if (j==k&&i~=j)

                    Hh(j+42,k+12,i+30)=V(i)*aY(i,j)*cos(Vth(i)-Vth(j)-Yth(i,j));      %thV

3 ~) ^5 t* n: S6 `3 p3 b8 M
elseif (j==k&&i==j)

                    temp=0;                   %thV

6 n5 t' d* P) Q1 J7 U. ffor l=1:30

                        temp=temp-aY(j,l)*V(l)*cos(Vth(j)-Vth(l)-Yth(j,l));   

; @0 y+ ]! i' y) hend

                    Hh(j+42,k+12,i+30)=temp+V(i)*aY(i,i)*cos(-Yth(i,i));

9 A$ c& D) X/ {' aelseif (j==i)

                    Hh(j+42,k+12,i+30)=-V(i)*aY(i,k)*cos(Vth(i)-Vth(k)-Yth(i,k));   %thV   

" o7 e5 k4 c4 S2 H% a  E$ a5 xelseif (k==i)

                    Hh(j+42,k+12,i+30)=V(j)*aY(i,j)*cos(Vth(i)-Vth(j)-Yth(i,j));  %thV

- j, u$ A$ V1 K& r" C$ \& z- }1 y2 send

, L  {8 z8 N9 B  x! b8 ^: pend

0 v+ {/ P9 S! r; Y, I$ [9 o1 x5 w* send

        Hh(13:42,43:72,i+30)=Hh(43:72,13:42,i+30)';

9 D# m: L, O$ W( D% P
end" s: V& k, ?/ h! Y  U3 B( e+ C
%ÖÁ´ËÒÑÐγɣ¨42-72£¬13-42£©ºÍ£¨13-42£¬43-72£©

. `; b1 S5 x# ?3 Y5 I7 W7 p6 L%HhÐγÉÍê±Ï

" j+ X9 k* F. \$ f+ h2 V2 H, M
%Fourth Step: Jg, Hg

    Jg=eye(42,42);

    Jg=[Jg;zeros(30,42)];

    Hg=zeros(72);

0 t) Z# Z1 _# G%Calculation Jacobian&Hessian matrix END

: A6 \* q/ M1 m
%%

; n9 M% k0 x! u$ u) t" }7 `6 ?
%Calculate Newton Iteration Îó²îµü´úÁ¿

# q% m7 w' m! \4 g5 P( K7 @
%Cal LX0-------------------------1

    LX0=Jf-Jh*Lam+Jg*(-MU_MIN+MU_MAX);

1 y, B% V7 K2 l+ w%Cal LLam-------------------------2

    LLam0=h;

    pferr=max(LLam0);

/ s  |' k  ~5 X: @5 h%Cal LMU_MIN-------------------------3

    LMU_MIN0=g-sl-gmin;

) j% I8 k" \, O# g$ l3 A. ^%Cal LMU_MAX-------------------------4

    LMU_MAX0=g+su-gmax;

# L0 t- X0 M' t  s9 ~3 l4 Z/ m+ m%Cal Lsl-------------------------5

    Lsl0=diag(MU_MIN)*diag(sl)*ones(length(sl),1)-MU*ones(length(sl),1);

$ o* w- d! G9 C' @* Q7 c% q1 M
%CAl Lsu-------------------------6

    Lsu0=diag(MU_MAX)*diag(su)*ones(length(su),1)-MU*ones(length(su),1);

* w2 O: a8 r7 N( I% @
%Calculate Newton Iteration Îó²îµü´úÁ¿ END!!!

) A' b  q) F0 y2 I8 L5 h1 Y%%

/ G+ t. w0 s$ J# }1 ]
%Calculate Newton Iteration ·ÂÉäÐÞÕýÁ¿

; u' }+ ~3 L. p( J$ _' G%1st Step: ·ÂÉäÐÞÕýÁ¿delXaf,·ÂÉäÐÞÕýÁ¿delLamaf

5 U( D! N( `+ F  V  m. \4 |%S1:H

     temp=0;  %֮ǰÒÑÓùýÁÙʱ±äÁ¿temp£¬ÔÚ´ËÇåÁã

" |1 [) v5 n" k$ h
for i=1:60

         temp=temp+Lam(i)*Hh(:,:,i);

0 `+ I7 w  i; T+ \6 I# Iend

     tempgMUg1=0;

     tempgMUg1=Jg*diag(MU_MIN)*Jg';

! a2 v* a; X4 j$ `
for i=1:42

         tempgMUg1(i,:)=tempgMUg1(i,:)/sl(i);        

6 u& ~0 i+ s' i  |# l8 [' {4 c$ s
end

     tempgMUg2=0;

     tempgMUg2=Jg*diag(MU_MAX)*Jg';

0 J$ C7 X" |% E! C
for i=1:42

             tempgMUg2(i,:)=tempgMUg2(i,:)/su(i);

2 ?  _5 o2 w$ r' A5 H5 }" f
end

     H=Hf-temp+tempgMUg1+tempgMUg2;

6 |+ I. a7 a$ Z* f
%S2:·ÂÉäKESAaf

     tempgMUg1=0;

     tempgMUg1=diag(MU_MIN)*sl+diag(MU_MIN)*LMU_MIN0;

4 T/ x1 T2 ]+ m- f+ tfor i=1:42

         tempgMUg1(i)=tempgMUg1(i)/sl(i);

1 I$ [+ A: A. t" T) F/ C; q
end

     tempgMUg1=Jg*tempgMUg1;

     tempgMUg2=0;

     tempgMUg2=diag(MU_MAX)*su-diag(MU_MAX)*LMU_MAX0;

+ v( g! N  Y, g8 y4 \% kfor i=1:42

             tempgMUg2(i)=tempgMUg2(i)/su(i);

; w; _3 x9 [' h. q& Z7 A
end

     tempgMUg2=Jg*tempgMUg2;

     KESAaf=LX0+tempgMUg1-tempgMUg2;

5 _* t5 Q0 S% ?" L  m. {
%S3:JACOBIAN

     JACOBIAN=[H -Jh;Jh' zeros(60)];

: L# o1 g! Z/ E3 u( a
%S4£ºRESULT

     RESULT=-[KESAaf;LLam0];

# z- K- u' D2 u1 b( T  p%S5:Cal ·ÂÉäÐÞÕýÁ¿delXaf,·ÅÉäÐÞÕýÁ¿delLamaf

     delXaf=-Jh'\LLam0;

     delLamaf=Jh\(H*delXaf+KESAaf);

%      temp=JACOBIAN\RESULT;

%      for i=1:72

%          delX(i)=temp(i);

%      end

%      for i=1:60

%          delLam(i)=temp(i+72);

%      end

' X1 D8 K! q& M! z& a) L5 c9 k

8 M* t6 g% X& ?
%2nd Step: ·ÅÉäËɳڱäÁ¿ÐÞÕýÁ¿delslaf, delsuaf

7 P$ o/ n. y# ^( c1 `3 N) z; [
%S1: delslaf

     delslaf=Jg'*delXaf+LMU_MIN0;

: T4 k9 D" e- I/ t' R%S2: delsuaf

     delsuaf=-Jg'*delXaf-LMU_MAX0;

* H5 `2 g: x) m6 a+ C

& i% G, b2 F/ }/ w' G0 }  p
%3rd Step: ·ÅÉäÀ ¸ñÀÊÈÕ³Ë×ÓÐÞÕýÁ¿delMU_MINaf,delMU_MAXaf

9 V- P" h! S7 o%S1: delMU_MINaf

     temp=0;

     temp=-Lsl0-diag(MU_MIN)*delslaf;

4 q) L# b7 Y7 V9 }/ o
for i=1:42

         temp(i)=temp(i)/sl(i);

( o: I& U' N8 O, d. V7 pend

     temp(13)=0;

     delMU_MINaf=temp;

. N& n7 N+ s9 B7 f: A2 p  g/ ^
%S2: delMU_MAXaf

     temp=0;

     temp=-Lsu0-diag(MU_MAX)*delsuaf;

# D% |) Z; s+ ^, ifor i=1:42

         temp(i)=temp(i)/su(i);

6 E5 @) y* _8 x1 a- ^5 wend

     temp(13)=0;

     delMU_MAXaf=temp;

: c* x+ z4 ?1 A5 o  @! G( n, X

& a% _! K$ P( s/ b* [%Calculate Newton Iteration ·ÂÉäÐÞÕýÁ¿ END!!!!!!!!!!!!!!

0 o  U4 c1 J0 z

2 Q5 p1 m% i4 ^5 _: b) @
%%

% \# s( \9 e% W4 [! d+ M. k
%¼ÆËã·ÂÉäÔ Ê¼²½³¤ºÍ¶Ôż²½³¤ STEPpaf£¬ STEPdaf

5 b. s* \/ |& p. a) K%S1: STEPpaf

4 f* c; h/ z1 G& y7 G
for i=1:42

+ H" w9 C: @/ B
if (delslaf(i)~=0)

            temp1=-sl(i)/delslaf(i);

' i" y2 S/ @+ ^5 u7 L" z" b
else

            temp1=Inf;

( {0 Y) G& H0 g" k" jend

7 G; K% t! x5 }, {% T6 z0 Uif (delsuaf(i)~=0)

            temp2=-su(i)/delsuaf(i);

# h8 Z/ d! Y! A0 D' C: u5 d. W3 }else

            temp2=Inf;

7 _+ K, }  w2 b5 j+ v- y, Q) Cend

* {4 f6 z' }; i& c6 Z. P8 ~3 |+ z9 {. j
if temp1<temp2

            min=temp1;

% |6 s" M8 E1 E, ?1 telse

            min=temp2;

+ O% h1 r1 }4 U, y: |' e' Y5 j! S
end

1 G  j6 T! {" \if min>1

            min=1;

  U7 K% y& }: R6 Z( P( Cend

0 `* B2 P" G% s0 ^4 d) K
end

    STEPpaf=0.9995*min;

5 a3 R3 y% f# \& f5 ]4 y$ O%S2: STEPdaf

: }# l  y" r2 D7 u& k, Z
for i=1:42

        temp1=-MU_MIN(i)/delMU_MINaf(i);

        temp2=-MU_MAX(i)/delMU_MAXaf(i);

4 J, D; \" t4 H, o* Wif temp1<temp2

            min=temp1;

( \/ T& R7 F3 y6 G! |else

            min=temp2;

" V1 Y9 l1 S4 F! X5 g5 Q
end

0 L, I& j! M3 B& w2 z6 x6 B8 Yif (i==13)

            min=0;

) u" Y. J4 G, l, C* }
end

% S$ N6 N9 a. R8 W
if min>1

            min=1;

; j+ d2 S/ m- m: Vend

  R# w  [2 ?. l5 k* N/ }3 N7 [end

    STEPdaf=0.9995*min;

2 Y- O: g8 f1 ^8 M%&frac14;&AElig;&Euml;&atilde;·&Acirc;&Eacute;&auml;&Ocirc; &Ecirc;&frac14;&sup2;&frac12;&sup3;¤&ordm;&Iacute;&para;&Ocirc;&Aring;&frac14;&sup2;&frac12;&sup3;¤ STEPp&pound;&not; STEPd END!!!!!

+ H) w2 j: j* v3 U

: Q8 {2 E" Z! B4 j8 L3 y%&frac14;&AElig;&Euml;&atilde;·&Acirc;&Eacute;&auml;&para;&Ocirc;&Aring;&frac14;&Ograve;ò&Ecirc;&yacute;&frac14;°&sup3;&Iacute;·&pound;&Ograve;ò×&Oacute;%

   ROUaf=(sl+STEPpaf*delslaf)'*(MU_MIN+STEPdaf*delMU_MINaf)+(su+STEPpaf*delsuaf)'*(MU_MAX+STEPdaf*delMU_MAXaf);

5 z/ ?+ {" M  h  J$ ~  I1 aif (ROUaf/ROU)^2<0.2

       SIGMA=(ROUaf/ROU)^2;

% k7 s; B5 P2 P! L$ t+ B- h
else

       SIGMA=0.2;

& v; T. G- D: u4 T2 \% l7 \
end

   MUaf=SIGMA*ROUaf/(2*length(sl));

. W" b8 l# Z6 _( _# j7 P%&frac14;&AElig;&Euml;&atilde;&Iacute;ê±&Iuml;%

   Lsl0=diag(MU_MIN)*diag(sl)*ones(length(sl),1)-MUaf*ones(length(sl),1)-diag(delslaf)*delMU_MINaf;

   Lsu0=diag(MU_MAX)*diag(su)*ones(length(su),1)-MUaf*ones(length(su),1)-diag(delsuaf)*delMU_MAXaf;

8 A$ K# c1 n" J) T, ]
%Calculate Newton Iteration &ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;

. G5 u. \3 Z2 n( P  L# K5 t& [
%1st Step: &ETH;&pound;&Otilde;&yacute;&ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;delX,&ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;delLam

+ l' U4 o3 L! i' s%S1:H

     temp=0;  %&Ouml;&reg;&Ccedil;°&Ograve;&Ntilde;&Oacute;&Atilde;&sup1;&yacute;&Aacute;&Ugrave;&Ecirc;±±&auml;&Aacute;&iquest;temp&pound;&not;&Ocirc;&Uacute;&acute;&Euml;&Ccedil;&aring;&Aacute;&atilde;

0 l$ H! N! \% ^( c+ H. a7 ]for i=1:60

         temp=temp+Lam(i)*Hh(:,:,i);

- h) y6 s  r) a1 B& D, ]  vend

     tempgMUg1=0;

     tempgMUg1=Jg*diag(MU_MIN)*Jg';

5 T5 F* B$ I# P- k
for i=1:42

         tempgMUg1(i,:)=tempgMUg1(i,:)/sl(i);        

: ~) C0 n6 O1 |+ e, B6 [: ]2 S3 |! V/ i
end

     tempgMUg2=0;

     tempgMUg2=Jg*diag(MU_MAX)*Jg';

0 K& [1 P+ N, k8 c+ Z5 g/ L, Pfor i=1:42

             tempgMUg2(i,:)=tempgMUg2(i,:)/su(i);

, a- q7 X  ]: n4 q4 |2 ^
end

     H=Hf-temp+tempgMUg1+tempgMUg2;

+ S& h. E/ Y  Z( h" s%S2:&ETH;&pound;&Otilde;&yacute;KESA

     tempgMUg1=0;

     tempgMUg1=Lsl0+diag(MU_MIN)*LMU_MIN0+diag(delslaf)*delMU_MINaf;

- O( \9 Q( J2 w/ hfor i=1:42

         tempgMUg1(i)=tempgMUg1(i)/sl(i);

0 ?7 ?$ u' a0 k& P5 _/ ?, ~end

     tempgMUg1=Jg*tempgMUg1;

     tempgMUg2=0;

     tempgMUg2=Lsu0-diag(MU_MAX)*LMU_MAX0+diag(delsuaf)*delMU_MAXaf;

7 o5 T7 S2 a5 g$ V# y7 ^0 Afor i=1:42

             tempgMUg2(i)=tempgMUg2(i)/su(i);

8 d# c% M  v3 t" R, Q4 x' bend

     tempgMUg2=Jg*tempgMUg2;

     KESA=LX0+tempgMUg1-tempgMUg2;

. z7 |& U0 s) O8 ?  d2 Y
%S3:JACOBIAN

     JACOBIAN=[H -Jh;Jh' zeros(60)];

, n( {1 g7 P& W% G8 U2 t& v
%S4&pound;&ordm;RESULT

     RESULT=-[KESA;LLam0];

6 V' ?, x- A, j7 D1 p1 [, _
%S5:Cal &ETH;&pound;&Otilde;&yacute;&ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;delX,&ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;delLam

     delX=-Jh'\LLam0;

     delLam=Jh\(H*delX+KESA);

%      temp=JACOBIAN\RESULT;

%      for i=1:72

%          delX(i)=temp(i);

%      end

%      for i=1:60

%          delLam(i)=temp(i+72);

%      end

9 o0 a9 Z- i' A' ~1 N

7 B9 j9 Q- \2 w" ?& W! R
%2nd Step: &ETH;&pound;&Otilde;&yacute;&Euml;&Eacute;&sup3;&Uacute;±&auml;&Aacute;&iquest;&ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;delsl, delsu

- u; o! w5 ]3 S$ i3 R
%S1: delsl

     delsl=Jg'*delX+LMU_MIN0;

5 w0 c# r4 w! B1 l. |4 W1 v%S2: delsu

     delsu=-Jg'*delX-LMU_MAX0;

/ K) n& C, K! ]0 n+ p2 [( E8 x( T  @% K

4 a) O) t6 L) `- `, F%3rd Step: &ETH;&pound;&Otilde;&yacute;&Agrave; &cedil;&ntilde;&Agrave;&Ecirc;&Egrave;&Otilde;&sup3;&Euml;×&Oacute;&ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;delMU_MIN,delMU_MAX

. [8 l' D; T* p5 }+ b! ?%S1: delMU_MIN

     temp=0;

     temp=-diag(MU_MIN)*sl+MUt*ones(42,1)-diag(delMU_MINaf)*delslaf-diag(MU_MIN)*delslaf;

  g2 T4 L# z0 a1 y
for i=1:42

         temp(i)=temp(i)/sl(i);

" N% C) ]  S3 x+ q
end

     temp(13)=0;

     delMU_MIN=temp;

* a1 D. M  T. f8 I) c%S2: delMU_MAX

     temp=0;

     temp=-diag(MU_MAX)*su+MUt*ones(42,1)-diag(delMU_MAXaf)*delsuaf-diag(MU_MAX)*delsuaf;

: e6 v# A, z# B) B# k" u  ofor i=1:42

         temp(i)=temp(i)/su(i);

) z9 g+ B9 l( m: r8 E0 F) z- [9 Pend

     temp(13)=0;

     delMU_MAX=temp;

2 h  l: \8 T) }) [7 f9 }/ r+ d

7 L; X+ d$ u% g% s
%Calculate Newton Iteration &ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest; END!!!!!!!!!!!!!!

7 a& o6 C' v, p%%

- F) l& u% d4 _
%&frac14;&AElig;&Euml;&atilde;&ETH;&pound;&Otilde;&yacute;&Ocirc; &Ecirc;&frac14;&sup2;&frac12;&sup3;¤&ordm;&Iacute;&para;&Ocirc;&Aring;&frac14;&sup2;&frac12;&sup3;¤ STEPp&pound;&not; STEPd

" S+ {5 S7 u( k: u
%S1: STEPp

/ r+ m; i. x0 S7 G6 z) e& afor i=1:42

6 K0 X) _& v: N9 j, Xif (delslaf(i)~=0)

            temp1=-sl(i)/delsl(i);

2 D- |8 R! g) p
else

            temp1=Inf;

7 ^6 N, ]% O4 `, \, e9 L1 i5 i
end

( x# N% N& v0 Y- E" l8 b, A. n2 Y
if (delsuaf(i)~=0)

            temp2=-su(i)/delsu(i);

  _& t/ d7 r  [else

            temp2=Inf;

9 K+ V6 i! w. Z9 O3 a
end

; ^1 g$ d( m* |! z  [! c; F0 rif temp1<temp2

            min=temp1;

# W0 B; V8 b6 R5 j1 Welse

            min=temp2;

7 j" Q* T  a9 \end

/ z! r% }9 T# P% N4 h' l8 ]if min>1

            min=1;

  F4 Z" ]0 x% T3 N$ v
end

5 W% u5 z7 S, x; c8 V! J4 Jend

    STEPp=0.9995*min;

3 Q/ g. D5 u! y0 z1 l+ r- I5 y- C  T
%S2: STEPd

9 D4 ~" c5 W( D* t, Lfor i=1:42

        temp1=-MU_MIN(i)/delMU_MIN(i);

        temp2=-MU_MAX(i)/delMU_MAX(i);

( m% |, r8 ^  p2 Q6 s# G( g, vif temp1<temp2

            min=temp1;

2 d+ A  n! z! R% Y( @. W/ J) F8 K- Relse

            min=temp2;

% H) Y$ J. |2 Zend

( x' G; w. w7 uif (i==13)

            min=0;

1 u, g( s: e: M2 ~5 R& H* oend

# z0 ^4 P2 Z7 G" B3 Tif min>1

            min=1;

3 P) d* O' b6 g4 s/ Hend

2 f1 L% y0 W* f) X' Yend

    STEPd=0.9995*min;

! a- X9 c/ Q5 a& K' f8 v& l3 r%&cedil;ü&ETH;&Acirc;&Ocirc; &Ecirc;&frac14;±&auml;&Aacute;&iquest;&ordm;&Iacute;&para;&Ocirc;&Aring;&frac14;±&auml;&Aacute;&iquest;

   X=X+STEPp*delX;

   Lam=Lam+STEPd*delLam;

   sl=sl+STEPp*delsl;

   su=su+STEPp*delsu;

   MU_MIN=MU_MIN+STEPd*delMU_MIN;

   MU_MAX=MU_MAX+STEPd*delMU_MAX;

   Pg=X(1:6);Qg=X(7:12);V=X(13:42);Vth=X(43:72);

- z; j+ i( S, W: R$ N+ h%&cedil;ü&ETH;&Acirc;&Ocirc; &Ecirc;&frac14;±&auml;&Aacute;&iquest;&ordm;&Iacute;&para;&Ocirc;&Aring;&frac14;±&auml;&Aacute;&iquest; END!!!

$ M) [: f* y# N0 b8 C& N0 x$ m

# p6 s2 B: ?, V
%&frac14;&AElig;&Euml;&atilde;&para;&Ocirc;&Aring;&frac14;&frac14;&auml;&Iuml;&para;%

   ROU=sl'*MU_MIN+su'*MU_MAX;

   MUt=SIGMA*ROU/(2*length(sl));

" U( r; t, ]: `0 P2 e% F6 r* l%%

1 R8 _8 i8 M- t- K%&Aring;&ETH;&para;&Iuml;&Ecirc;&Ccedil;·&ntilde;&sup3;&not;&sup1;&yacute;&micro;ü&acute;ú&acute;&Icirc;&Ecirc;&yacute;

    ik=ik+1;

7 S6 @, o+ o: ^/ L% p0 K
if ik>IterNumMAX

        disp('IterNum ERROR!!!!!');

# m0 I( U) c) I1 j* }
break;

) W# W$ `2 @# |8 U$ b
end

* I- I* n/ E& I. M( t

end

et = etime(clock, t1);

%%

%&Ecirc;&auml;&sup3;&ouml;&sup2;&iquest;·&Ouml;

* p; r4 @; h; Q8 _' B7 x

if (ik<=IterNumMAX)

$ N- u; ?1 j! c0 X, ?8 V%%

3 f* P3 [  R8 V( e" ?3 D$ b8 z%&Ccedil;ó&micro;&Atilde;×&icirc;&acute;ó×&icirc;&ETH;&iexcl;&micro;&Egrave;&Ecirc;&frac12;&Icirc;ó&sup2;&icirc;

    F=0;

for i=1:6

    F=F+gencost(i,5)*(X(i)*baseMVA)^2+gencost(i,6)*X(i)*baseMVA;

end

8 P% k! b( e* C) Z# `

max=-10;min=10;

for i=1:ik

    temp=LLam0(i);

) c# D/ B' w4 R6 c5 n. pif (temp>max) max=temp;

; ?) _) \6 `- D/ V5 l" Cend

) T+ j  I9 m) N0 _- V; @$ Bif (temp<min) min=temp;

. _4 |5 K" L) b: `
end

end

max

ik

%%

%&Ecirc;&auml;&sup3;&ouml;&frac12;á&sup1;&ucirc;

    fprintf('\nConverged in %.2f seconds', et);

    fprintf('\nObjective Function Value = %.2f $/hr', F);

    fprintf('\n================================================================================');

    fprintf('\n|     Bus Data                                                                 |');

    fprintf('\n================================================================================');

    fprintf('\n Bus      Voltage          Generation             Load        ');

    fprintf('  Lambda($/MVA-hr)');

    fprintf('\n  #   Mag(pu) Ang(deg)   P (MW)   Q (MVAr)   P (MW)   Q (MVAr)');

    fprintf('     P        Q   ');

    fprintf('\n----- ------- --------  --------  --------  --------  --------');

    fprintf('  -------  -------');

+ v. _0 c( W6 w) a
for i = 1:30

        fprintf('\n%5d%7.3f%9.3f', i,X(i+12),X(i+42)*180/pi);

) n7 R. p  D$ s4 m) A
if (i<=6)

            fprintf('%10.2f%10.2f', X(i)*baseMVA, X(i+6)*baseMVA);

  N: d6 v2 `$ d( Qelse

            fprintf('       -         -  ');

6 X) S( L; R/ H! d
end

        fprintf('%10.2f %9.2f', bus(i, 3) , bus(i, 4) );

        fprintf('%9.3f %9.3f', Lam(i),Lam(i+30));

' E0 l7 o) R; ]. Q, zif (i==30)

            fprintf('\n                        --------  --------  --------  --------');

            fprintf('\n               Total: %9.2f %9.2f',(X(1)+X(2)+X(3)+X(4)+X(5)+X(6))*baseMVA,(X(7)+X(8)+X(9)+X(10)+X(11)+X(12))*baseMVA)

            fprintf('\n');

$ B+ X( ?' i+ Iend

. l% w/ S( @. c& F" ^end

! c) c/ W% m% p* f%%

%%&Ecirc;&auml;&sup3;&ouml;&Ccedil;ú&Iuml;&szlig;&Auml;&acirc;&ordm;&Iuml;main

iterNum=[1:ik]';

%axis([0,ik+1,min,max]);

hold on;

f = fit(iterNum,errArr,'spline');

f=feval(f,iterNum);

plot(iterNum,errArr,'o',iterNum,f,'-');

end