[讨论] 关于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)

  l' O5 U& a' h. B0 u& L3 }%%

* A* B% }7 I) @" b6 |' P
%Calcute h,g matrix

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

        errArr=[errArr;ROU;];

    SIGMA=0;

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

$ U7 G0 c- a4 B# n) nfor i=1:30

        temp=0;

2 }! {/ V7 K8 s& nfor j=1:30

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

  i  q8 {0 b1 ~& q. Q' h3 |
end

4 U9 E& Q) Y  Y" |9 T, h
if (i>6)

            tPg=0;

$ C- u0 i4 I! h8 [else

            tPg=Pg(i);

) B) }; W: t$ Z- B5 lend

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

8 @2 e3 E7 A4 {
end

+ M( l6 d- ?" V  W: T
for i=1:30

        temp=0;

) a$ Q5 W% |6 t; M5 L. X/ i; @% [- J
for j=1:30

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

8 ~( ]* u& z# P3 r. l6 e/ hend

; u" V, A* O! {$ k( Mif (i>6)

            tQg=0;

* m6 [+ Z: n& g) o" m2 l; l, o8 B
else

            tQg=Qg(i);

$ S$ O* t0 N5 M& o/ C9 g- }end

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

# {; R5 G/ x! u% l7 I. mend/ W2 l4 l2 A4 j: {  U- S. E; d
% Cal  h END

6 j/ A7 O8 z' v, S! v# ?* f

7 V' n' Q# d7 v1 J. Y* V' n; c. v+ x6 k
for i=1:6

        g(i)=Pg(i);

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

% B/ k3 L+ l/ r4 _% v8 o4 Uend

% p( D; Z+ H  R5 e- q: J" G' Z
for i=1:30

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

2 Y( }5 t/ u) x( m) Q' _' {( mend+ r# w* z2 l7 Z/ Z* [+ l
% Cal  g  END

6 Q+ }7 z" {' x& H% Q$ D0 P0 d
%Calcute h,g matrix END

1 T  Z( l! [! V9 F%%

5 P8 g. a; z  W* p! \! O. t7 a%Calculate Jacobian&Hessian matix

5 j' D( Z+ s) Z# i, h" y%First Step: Jf,Hf

# m# \- Z5 W& S
for i=1:6

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

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

% I- v8 `) n+ S* g2 x& A1 J0 eend

  i! L% Q6 f! I/ X- Y& L! C; \) ?
%Second Step: Jh, hΪµÈʽԼÊø

2 {/ D& W# H6 Q& X0 t# I2 B
for i=1:6         %Ç°6ÐжÔPgÇóµ¼£¬ÓÉ´ËÒÑÇó³ö

        Jh(i,i)=1;

6 w5 D9 d4 N# B
end

( R: h. U; f% F* k4 R. y3 u+ |for i=7:12        %7-12ÐжÔQgÇóµ¼£¬ÓÉ´ËÒÑÇó³ö

        Jh(i,i+24)=1;

8 J4 ^' C5 W# A! @end

4 R. a) E) J1 _+ }for i=1:30       %ÐγÉ13-42ÐеÄ1-60ÁÐ

& t( Y4 N; R& k  C- S1 b1 mfor j=1:30

            tempVp=0;

            tempVq=0;

6 \: [& V+ G% z* d7 qif (j==i)

% P; F% [7 Q: }. |, Pfor 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));

2 _0 U/ C; c7 e& e; M, ?0 e# C- i( z
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));

) a9 }' k# x2 i5 T+ u/ ^else

                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));

" b, u- S5 p: R! w' P, }end

$ S0 G- r  d- |. ]4 c$ @( Tend

0 u0 D) o6 s) [! c, p) r
end

8 K/ b( r. S7 Z3 ~
for i=1:30       %ÐγÉ43-72ÐеÄ1-60ÁÐ

* k0 I# R9 u6 z8 ^: }for j=1:30

            tempVp=0;

            tempVq=0;

, o( A- l8 v1 @" U# Y; f% B8 Gif (j==i)

5 r0 j: ~* U; P  b! G
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));

' n7 g9 r. o$ s; E
end

                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;

* g% Q5 J6 B, a4 c6 q2 qelse

                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));

5 F5 \" ~! ^; Q4 {
end

0 t2 d" h6 G9 Y' u. ~0 A3 zend

2 C( ?7 C, j; M( P& _
end

5 k+ h4 p# K  d  a
%Third Step: Hh

8 X. p1 R% @3 }3 a  b6 D%Óй¦²¿·Ö

/ o1 \& i, i9 I# R2 i0 g, y
for i=1:30

8 [0 j- B8 B0 l- V- Y: u( F
for j=1:30

* [& N3 m& r5 Lfor k=j:30

1 |, v, t( o5 @3 b
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

. E0 w/ R% E9 B/ v% j% V! |: Felseif (j==k&&i==j)

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

                    temp=0;                   %thth

$ K. v* _- z2 U7 S
for l=1:30

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

! h& f" O6 F$ |2 V9 iend

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

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

( x5 I8 Z( C( `2 k- v5 B) Q
elseif (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);

$ n- W. `& G3 Celseif (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);

5 V* `/ P2 [9 U+ `end

2 {' T: F9 u  h5 n6 w: A
end

$ [+ p, i+ X$ F3 I/ j; R
end

" H2 C4 W% }9 Q, M
end
1 B& c) U" O. _( F. S7 `1 Z%ÖÁ´ËÒÑÐγɣ¨13-42£¬13-42£©ºÍ£¨42-72£¬43-72£©

5 V" s# x% D6 z- ?0 i
for i=1:30

  z' u% q+ N& ?/ }$ _
for j=1:30

  d; a) |5 T& `# Qfor k=1:30

) \4 g; o9 f2 Z% q
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

; [$ b' H% M  l& k% R8 i, h# w( a
elseif (j==k&&i==j)

                    temp=0;                   %thV

+ k, y0 H' F2 m+ P" j! [% v$ v
for l=1:30

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

% u( Q1 E- r* m4 ~# V/ nend

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

& V5 n2 ~' @& d9 F9 r  ?9 pelseif (j==i)

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

; l( u7 f2 H, p7 c# S7 I
elseif (k==i)

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

% @. N4 R1 R& q% u
end

* h9 i/ R. ?+ W, @1 Q* Mend

* s- _9 {, m3 D4 c& l% [4 S" aend

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

, ]) C0 A4 e6 c2 L9 y8 Jend
; c1 a8 b4 ]1 \' L& m/ [* l3 ^3 ?%ÖÁ´ËÒÑÐγɣ¨42-72£¬13-42£©ºÍ£¨13-42£¬43-72£©

+ ~" y* A7 ]" D+ A+ L
%ÎÞ¹¦²¿·Ö

' u9 g# `- Q% N$ T) q+ l& ^
for i=1:30

3 {" |% Y6 w+ qfor j=1:30

# N/ |2 X; p. `( X2 vfor k=j:30

( Y9 P  c# R) z7 R6 n
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

* U# X. d9 O( ?! Oelseif (j==k&&i==j)

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

                    temp=0;                   %thth

% g7 R5 t$ H& ?6 b$ F  Q6 |for l=1:30

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

& ?+ \. u. h( P2 i
end

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

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

1 U% G3 S% J1 m! t- v8 `
elseif (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);

' B. p7 v  l1 i. H" B0 `8 d" |4 z0 m
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);

2 E6 c8 z) X" u; g# b) [end

) N0 \6 J. r/ ?end

8 c, Q: a/ n1 D6 x2 Y, cend

; E% \1 l# B2 _5 f% D! \2 G
end+ I8 R( G2 h  Q. }% D
%ÖÁ´ËÒÑÐγɣ¨13-42£¬13-42£©ºÍ£¨42-72£¬43-72£©

" k1 O6 U  h1 W2 h) k& r3 X
for i=1:30

" D4 i. B, B; b; Sfor j=1:30

4 E1 Y- X: i5 z& z' e$ [
for k=1:30

) q& P- \+ B0 ?# q
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

& u* }) h( g7 X0 l* N2 H. o0 delseif (j==k&&i==j)

                    temp=0;                   %thV

  W& V4 w1 n% Y% n+ ]% t) m% t
for l=1:30

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

. B1 T( w9 z3 _: w* Z3 J' m
end

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

! m- ^" W% m* Yelseif (j==i)

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

! Q1 P! a! h. x2 R: L6 k
elseif (k==i)

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

6 p) E& `5 Q3 R! w& p
end

3 J! v* K; U. Z3 ~6 Zend

0 e% w/ N9 D. send

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

" B/ Q2 W+ `6 T) k; f$ M2 l( C/ yend# F" c7 Q" C9 ^" ]/ O0 D
%ÖÁ´ËÒÑÐγɣ¨42-72£¬13-42£©ºÍ£¨13-42£¬43-72£©

& u  o* z* E0 E; ?- H6 c/ g$ \
%HhÐγÉÍê±Ï

. T4 _' F, I" Y, n9 a, P
%Fourth Step: Jg, Hg

    Jg=eye(42,42);

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

    Hg=zeros(72);

0 p. e" j2 ]0 d, T6 G) F' h
%Calculation Jacobian&Hessian matrix END

# c4 z6 c) m- y* _* O%%

( G2 K% j; L0 |/ H$ O. o
%Calculate Newton Iteration Îó²îµü´úÁ¿

# `8 x3 N7 b, x, E
%Cal LX0-------------------------1

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

9 g6 J+ h+ y  }9 [%Cal LLam-------------------------2

    LLam0=h;

    pferr=max(LLam0);

4 I6 m( I2 e9 s8 ?/ [%Cal LMU_MIN-------------------------3

    LMU_MIN0=g-sl-gmin;

  ]) U( t- g2 \* m' H%Cal LMU_MAX-------------------------4

    LMU_MAX0=g+su-gmax;

+ v/ ~$ s0 H" L
%Cal Lsl-------------------------5

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

# U4 o! C) K. ~. S%CAl Lsu-------------------------6

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

9 \; y) \( q8 s8 `7 b%Calculate Newton Iteration Îó²îµü´úÁ¿ END!!!

9 u; d+ [3 T/ Y: T  T4 c- x1 A  [
%%

, s% s4 }& I- p( x& |8 M% [2 L) d
%Calculate Newton Iteration ·ÂÉäÐÞÕýÁ¿

( c8 o/ g0 W( P3 T' Q%1st Step: ·ÂÉäÐÞÕýÁ¿delXaf,·ÂÉäÐÞÕýÁ¿delLamaf

' V: z! e2 e) D$ p
%S1:H

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

, J7 U6 ~* G! q  A0 U  {for i=1:60

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

/ @5 O! O3 L  ~% P. A3 c0 l- Aend

     tempgMUg1=0;

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

; a, C6 O$ b5 m
for i=1:42

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

4 J9 Y6 i2 _2 p' ]7 c" f2 s, j9 send

     tempgMUg2=0;

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

. F& I  E$ f7 F, J
for i=1:42

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

6 d& r( N) l1 Y6 Oend

     H=Hf-temp+tempgMUg1+tempgMUg2;

# F$ J3 {: a, q
%S2:·ÂÉäKESAaf

     tempgMUg1=0;

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

% _" X. Q9 k9 x' l0 @for i=1:42

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

+ l5 o1 ^) p! J' p6 R" z$ X' ]
end

     tempgMUg1=Jg*tempgMUg1;

     tempgMUg2=0;

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

% n. T1 }) s) ?; z  U# ]9 i
for i=1:42

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

6 |5 e& u3 f# T; m8 `! D/ a- jend

     tempgMUg2=Jg*tempgMUg2;

     KESAaf=LX0+tempgMUg1-tempgMUg2;

. t2 l& D% E$ K. a1 V7 f: a9 Q3 f%S3:JACOBIAN

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

' W* r% G- X5 n3 f- v%S4£ºRESULT

     RESULT=-[KESAaf;LLam0];

! l9 V: i$ n1 D$ `
%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

+ @, i  J/ l  q! a1 u

5 X# M# l4 c8 d* d* n- ]3 }6 W%2nd Step: ·ÅÉäËɳڱäÁ¿ÐÞÕýÁ¿delslaf, delsuaf

2 F2 b8 h6 r: }: V/ S%S1: delslaf

     delslaf=Jg'*delXaf+LMU_MIN0;

8 [+ o8 G( T9 U7 ~, \%S2: delsuaf

     delsuaf=-Jg'*delXaf-LMU_MAX0;

5 E/ E6 p) |! g/ O% r3 Z2 n5 X+ j- j

- Y# F) c; J! r- ?( i/ ?4 j7 R+ ^. S5 k%3rd Step: ·ÅÉäÀ ¸ñÀÊÈÕ³Ë×ÓÐÞÕýÁ¿delMU_MINaf,delMU_MAXaf

$ L' @2 F& Q% V%S1: delMU_MINaf

     temp=0;

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

5 F+ R6 J% R( J9 I
for i=1:42

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

, [# G+ k( O* Z+ `' O! i0 ~
end

     temp(13)=0;

     delMU_MINaf=temp;

8 j: x$ |& b& w2 |* y/ D%S2: delMU_MAXaf

     temp=0;

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

4 h9 U1 d) Z% f! ]( _( o& Q
for i=1:42

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

) ~7 I7 R5 v: Q$ ~; l: @end

     temp(13)=0;

     delMU_MAXaf=temp;

' p6 o" a4 G/ }; B9 y3 i

/ s( a9 a% G! |# Y" M%Calculate Newton Iteration ·ÂÉäÐÞÕýÁ¿ END!!!!!!!!!!!!!!

: j( w5 P3 H) u0 l, x

! Y. m2 u* |' d5 ?% E6 k%%

# x0 K5 n# @( E* `: ]
%¼ÆËã·ÂÉäÔ Ê¼²½³¤ºÍ¶Ôż²½³¤ STEPpaf£¬ STEPdaf

3 [$ G, O) {% Q6 x%S1: STEPpaf

% f" Z$ q# S# k; ~
for i=1:42

( u. a+ j. k- p: A0 M$ {0 G5 n( oif (delslaf(i)~=0)

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

/ ]# x' }1 {7 G! s
else

            temp1=Inf;

/ M0 a% ^7 l5 R4 c% i8 n
end

/ Z* l5 F/ C* ~4 j' @if (delsuaf(i)~=0)

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

. S9 m+ r2 R* Q. g6 T6 Y
else

            temp2=Inf;

. J8 K6 ~! C9 O' ?6 ?
end

! B! ]9 d$ ^! q; @: w8 cif temp1<temp2

            min=temp1;

) U6 }7 P/ ^$ |
else

            min=temp2;

" Y% I8 `7 H' q$ \  h) Q& Z
end

, T( l4 f0 @: m9 e4 J1 s
if min>1

            min=1;

$ \1 C- i$ y/ J' Bend

, y0 C. o' h# ?; v" O, h8 Oend

    STEPpaf=0.9995*min;

6 d0 p3 [" k) \# N) m%S2: STEPdaf

- @: N3 e' V- B3 k; t, o7 T* N
for i=1:42

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

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

7 w( v9 e6 B5 x- X
if temp1<temp2

            min=temp1;

% w1 Y$ Q4 T- G$ ~* S
else

            min=temp2;

- o" W) Q2 Q/ g( G& c; bend

5 d5 |3 @# g7 i" a. S
if (i==13)

            min=0;

+ y* C8 N5 i; t, ~8 H4 tend

- T2 u2 L# C; }, N' C9 W9 Z- Mif min>1

            min=1;

& N6 j' Y3 W$ {2 Jend

2 L+ q& g+ Y) M# v" c8 Gend

    STEPdaf=0.9995*min;

! {( |8 U5 V. d, |4 [2 t7 X3 [%&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!!!!!

0 Y# R  m4 n( U- K# h, F, F

8 g1 _  E0 H4 ~( q%&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);

# L/ M# b0 B* |3 |4 h
if (ROUaf/ROU)^2<0.2

       SIGMA=(ROUaf/ROU)^2;

% r: B& B. T  V3 Z- A3 b# ?else

       SIGMA=0.2;

- n9 @8 u1 }/ s
end

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

$ K% T0 g$ D5 Z8 h! P6 Z%&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;

/ {% k' V  Z# ?%Calculate Newton Iteration &ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;

3 t9 K% q# x9 ^/ x+ E) ?; [6 }+ r) i%1st Step: &ETH;&pound;&Otilde;&yacute;&ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;delX,&ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;delLam

( b7 S6 K9 n7 U' t! Q1 j
%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;

' L, T+ ^* I4 U# b/ j) ~for i=1:60

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

! {0 |" ^, d' `9 |# i3 N& i' |
end

     tempgMUg1=0;

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

' W: ?: D" w9 A5 A! z4 Ifor i=1:42

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

  j3 E) j8 y) u4 J8 \/ `end

     tempgMUg2=0;

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

& ~! D- e" k9 |; M
for i=1:42

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

5 X" l7 `; W3 ]! Q& h1 z' K% H, x5 Rend

     H=Hf-temp+tempgMUg1+tempgMUg2;

( l! w5 ]2 D5 G% ]6 a%S2:&ETH;&pound;&Otilde;&yacute;KESA

     tempgMUg1=0;

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

- ~8 l( ~4 d/ z& [
for i=1:42

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

! d( ~" Z/ b6 [: g8 Z2 r! hend

     tempgMUg1=Jg*tempgMUg1;

     tempgMUg2=0;

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

9 S  S. i% v$ y4 ~
for i=1:42

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

1 H1 S  o. l' y) E, v
end

     tempgMUg2=Jg*tempgMUg2;

     KESA=LX0+tempgMUg1-tempgMUg2;

+ h5 W( L# y& _' J2 L" G4 N5 T+ [
%S3:JACOBIAN

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

+ h7 c/ ^  `: x' C8 I
%S4&pound;&ordm;RESULT

     RESULT=-[KESA;LLam0];

. D5 e/ D& x- |) `2 N6 w
%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

8 w" B# O8 b9 D6 l' C, ]" A  H

3 x1 Z% i; L* a  p/ ~
%2nd Step: &ETH;&pound;&Otilde;&yacute;&Euml;&Eacute;&sup3;&Uacute;±&auml;&Aacute;&iquest;&ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;delsl, delsu

- @0 T0 a% a! W% S: b" A. b0 R
%S1: delsl

     delsl=Jg'*delX+LMU_MIN0;

/ C+ l3 ?+ h) V5 v7 N' i' u9 M) j%S2: delsu

     delsu=-Jg'*delX-LMU_MAX0;

( F" P* c  U& ^' V3 y

7 l$ k0 E% a1 p' W* X! |* [. c
%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

% w7 r% d+ m+ v, B' N5 r& Z4 c
%S1: delMU_MIN

     temp=0;

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

6 h1 R" D! h- n$ Ofor i=1:42

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

6 s3 y% r) V, j, p1 a5 |+ _
end

     temp(13)=0;

     delMU_MIN=temp;

  B, {8 G0 {/ j
%S2: delMU_MAX

     temp=0;

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

  q# R, `1 }4 Q. y
for i=1:42

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

3 f8 Z0 {" U. D( T' _end

     temp(13)=0;

     delMU_MAX=temp;

- V# e$ Y1 L" h$ V9 O

) ^$ U7 V, V: o7 s5 \
%Calculate Newton Iteration &ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest; END!!!!!!!!!!!!!!

7 y8 ?1 V' u. B/ K; H- O) M
%%

; h! H+ h; p3 v) ?2 m, n5 L4 r4 W$ H%&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

! ^7 Y% t3 @% I; ?5 b
%S1: STEPp

$ v- j+ {' m  g2 v1 E
for i=1:42

' v# i+ v) Q" J, Q( c5 n
if (delslaf(i)~=0)

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

8 V0 n7 {: ~$ l5 e9 h
else

            temp1=Inf;

/ O: @7 q/ P: c/ a. p  B! |end

3 e6 a  M" T) k  Y2 y# O  rif (delsuaf(i)~=0)

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

4 E% u5 w; l9 l9 ]" Uelse

            temp2=Inf;

2 R8 c# ?+ }" @
end

+ f( t0 _! r2 r+ i% f7 V. ]) eif temp1<temp2

            min=temp1;

9 }2 h4 `8 a: y, r7 welse

            min=temp2;

5 ^' ?+ p1 `. b+ e' M5 Q+ s9 G6 qend

" d2 K7 }. \  e" }* r; \- M
if min>1

            min=1;

+ z/ r0 |# w' n) d" w, Z! nend

3 C" F! S/ y& a- n( D) A% X' n4 g
end

    STEPp=0.9995*min;

$ I: T+ B9 l6 G; f%S2: STEPd

, N8 m, S- }( g
for i=1:42

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

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

) t2 l, P. g2 {+ a6 R" m2 k
if temp1<temp2

            min=temp1;

, v* J, X6 m4 g' d
else

            min=temp2;

2 F: B9 D# M$ @* m. O- g* C
end

: W1 O" `8 [$ B. A2 c- J
if (i==13)

            min=0;

( z$ X. ]" Z: t" ?$ yend

$ d+ \5 ~& P2 O/ u
if min>1

            min=1;

. e, p% U- i7 w2 d* f, I7 B7 x$ y
end

- N  S8 Z3 V$ ]* u
end

    STEPd=0.9995*min;

' S" h) z# n/ a; u
%&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);

* E% c/ [2 o& C1 D%&cedil;ü&ETH;&Acirc;&Ocirc; &Ecirc;&frac14;±&auml;&Aacute;&iquest;&ordm;&Iacute;&para;&Ocirc;&Aring;&frac14;±&auml;&Aacute;&iquest; END!!!

9 w2 D- B$ @: y. R) I9 @

* z5 X% g; _+ t+ E( f) t) P+ R: W. {- H
%&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));

- w. \, m7 E6 t  _: _+ h%%

0 j0 E% C; D1 J' x
%&Aring;&ETH;&para;&Iuml;&Ecirc;&Ccedil;·&ntilde;&sup3;&not;&sup1;&yacute;&micro;ü&acute;ú&acute;&Icirc;&Ecirc;&yacute;

    ik=ik+1;

# w9 x; \/ u% ?: e; x6 b
if ik>IterNumMAX

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

9 V. l% l( A- m
break;

' f' _5 V) r7 }" ~5 Iend

: d7 P/ {* {6 Y9 c

end

et = etime(clock, t1);

%%

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

3 Q9 n7 \/ N" E

if (ik<=IterNumMAX)

$ Q6 ]# i# J6 n+ `%%

9 X6 N2 t0 P" h
%&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

' |7 e9 Y: `& O. i2 [

max=-10;min=10;

for i=1:ik

    temp=LLam0(i);

' r9 }, K- _6 ]0 [. @
if (temp>max) max=temp;

9 _+ X( N/ J8 c+ @- H- z, ^& ~0 b
end

3 t0 m5 b1 n4 \  Y6 S' `+ q# |if (temp<min) min=temp;

; e) q; j8 V/ p+ t& S% r) W# g, Gend

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('  -------  -------');

/ E# T$ k8 ~" S6 I. M/ Wfor i = 1:30

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

; W9 I0 B- j4 p0 ]8 H/ E5 ~
if (i<=6)

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

% Q! D5 |" ?- ]! l& |. u  ~
else

            fprintf('       -         -  ');

5 u4 r5 H( E' u9 r7 }end

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

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

( c! i8 J; c( y# J0 ~( @
if (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');

2 X7 X" c8 d! M% E. Q- Q! wend

3 O0 g7 Y" E) t7 o7 S& Uend

' a, p. `9 ]1 C- t/ h  ^" M' Z%%

%%&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