关于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)

) E" z; x* Q8 \- l
%%

5 P1 ~/ A, a0 ]4 C7 |/ H%Calcute h,g matrix

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

        errArr=[errArr;ROU;];

    SIGMA=0;

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

# j2 J1 O' J4 r* b3 x* }( Tfor i=1:30

        temp=0;

, d' R8 V3 J2 U* }( U& `for j=1:30

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

5 |( h$ {0 G$ ~3 q$ J% E
end

" b2 p* S1 I. `1 E0 ]+ j# a4 |' ^if (i>6)

            tPg=0;

  C5 Y% p0 N8 g0 K1 T5 ^! l) o
else

            tPg=Pg(i);

% D3 e$ M7 h: H
end

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

- O& O* x2 s8 V: c# P% Y
end

8 L. ^* c0 Z  h, S' @$ K: Tfor i=1:30

        temp=0;

0 V9 \4 @" v: X; W0 x0 dfor j=1:30

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

2 i6 I3 _1 e0 U2 f- S' e/ u$ O' Rend

% j7 ]; C5 i; _: S( Yif (i>6)

            tQg=0;

$ v2 O9 T0 z8 F, R, T; Z  x. r
else

            tQg=Qg(i);

! ~; z& q' g6 V, M" a2 S% p$ C
end

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

" P# K5 x) q: }+ k7 D8 W* }end
; N( R$ |+ r+ o4 y- I6 ], F% Cal  h END

7 u# x4 C5 w+ t$ f3 z' X' Y

6 s) j+ ?3 d6 E9 O9 R5 Dfor i=1:6

        g(i)=Pg(i);

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

8 _2 k  t) J/ c- O1 c
end

' b3 w( l, ]7 @" A
for i=1:30

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

0 F0 G3 S* L3 N
end
* g. P/ |1 w6 I% Cal  g  END

' k3 ]  B& i4 K# d
%Calcute h,g matrix END

  A6 c3 M( t! u0 r: f" B
%%

4 U* _1 U) n) m% V) i* [7 j
%Calculate Jacobian&Hessian matix

9 n* p' y. J7 R, T) q7 ^4 e6 \. @
%First Step: Jf,Hf

  g1 g5 @2 H( d+ C2 D! _
for i=1:6

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

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

. j3 V( E; l, L4 U" Q  Q3 T/ R5 oend

' F1 b+ P- ^. Z7 r' J) _# S%Second Step: Jh, hΪµÈʽԼÊø

4 s8 B& [: l7 ^/ ~7 a
for i=1:6         %Ç°6ÐжÔPgÇóµ¼£¬ÓÉ´ËÒÑÇó³ö

        Jh(i,i)=1;

, u6 I% A4 _( I5 dend

- R' d/ k4 U0 ], z5 u& n+ b( p9 Qfor i=7:12        %7-12ÐжÔQgÇóµ¼£¬ÓÉ´ËÒÑÇó³ö

        Jh(i,i+24)=1;

+ O. i8 e1 D9 q/ N& C9 k& V$ jend

- r/ J) F9 Q% k; o
for i=1:30       %ÐγÉ13-42ÐеÄ1-60ÁÐ

9 Q! L  [1 f( O9 q3 ifor j=1:30

            tempVp=0;

            tempVq=0;

( S+ e' ~) y* S9 q- H4 Cif (j==i)

% U6 u/ r( k3 w# k: y0 Ffor 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));

5 n' Q0 {* Y% G% \' _7 f
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));

1 x" s7 a8 r; }* p/ O" helse

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

2 `5 u. E0 v1 M0 Q2 o- T' V
end

* c" p$ A- K. }9 {! j$ U7 T& \
end

5 W% p) u0 [7 l4 k& u( r
end

; T& x5 S$ c' d
for i=1:30       %ÐγÉ43-72ÐеÄ1-60ÁÐ

) U# f3 ~6 w+ V: [1 U
for j=1:30

            tempVp=0;

            tempVq=0;

0 O1 V7 t% s; ~  |5 ]2 _1 L- f
if (j==i)

/ K9 Z6 ~1 t4 ?0 ~2 Q, e5 @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));

8 [$ p! X4 ?! U9 t9 Z
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;

4 z* U) _# z2 C5 i( U* j; C6 H
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));

) z; X2 k! q4 T* Hend

1 L5 {' A) z) {8 E0 P" k3 P( M
end

! s; w3 f9 a5 X" P: C( pend

- |* ]5 z: C/ b
%Third Step: Hh

* u) a3 D" H( f8 n%Óй¦²¿·Ö

1 \7 j. ^6 q9 |. u3 L8 qfor i=1:30

) w1 c! W! Q3 r8 T6 Mfor j=1:30

6 c) X& @: E1 S# G, k* u
for k=j:30

- w( L$ W* f" {0 i$ O5 h# i
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

+ }$ Q: M- R- m1 L$ T
elseif (j==k&&i==j)

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

                    temp=0;                   %thth

+ R: e6 {, a$ Y" c% ~6 p% @4 Q
for l=1:30

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

: w/ |' U- a& l- ~* w, wend

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

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

0 m! p2 ]$ f; o' oelseif (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);

" I+ u6 w5 N  M
elseif (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);

8 ^& o4 P4 T/ S) [end

+ l$ P7 w" e. R/ H7 }& w/ G
end

" Y/ @5 {( @1 b2 _& R3 b2 T
end

7 r% }5 H' G0 F& tend
& A0 A" y3 _1 J$ O7 t6 o2 T9 G%ÖÁ´ËÒÑÐγɣ¨13-42£¬13-42£©ºÍ£¨42-72£¬43-72£©

# \2 `- a. n( D) u# P( X9 y. gfor i=1:30

7 y' r! j% Z9 F" o+ x2 b& z
for j=1:30

# g8 j6 r# s  e: Y' Lfor k=1:30

2 S6 Y$ f* {3 k" y4 k5 Qif (j==k&&i~=j)

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

9 i; X8 A1 R; N% G3 ?- E& @8 f4 \
elseif (j==k&&i==j)

                    temp=0;                   %thV

8 _& W: P/ S! B& f9 R* T* T
for l=1:30

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

) X$ K7 y! b. I# r4 q9 N5 B* cend

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

: J$ {; @( ~5 v% I$ w' t4 z- Q1 nelseif (j==i)

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

* l! k  s0 S/ P$ W6 A7 F5 k
elseif (k==i)

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

2 C& T0 G9 t% tend

+ y) i" k) k* c/ f3 B2 p3 g
end

, R# j8 S& D/ M
end

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

* @0 g7 d6 W4 i: R7 ?2 p
end4 A, j4 U) w6 E0 ?+ d# Y% Q
%ÖÁ´ËÒÑÐγɣ¨42-72£¬13-42£©ºÍ£¨13-42£¬43-72£©

4 T; h7 h, H) m+ e%ÎÞ¹¦²¿·Ö

6 ?" c$ C( c- N' [! {' Mfor i=1:30

) E8 v; Q. _4 Z( n& A6 D
for j=1:30

: c  M* j& g* b( J( W, I
for k=j:30

8 j* Q6 O3 B  ?0 ^" W7 r2 Dif (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

( n* J' i  i9 C( W( N7 h/ Y, Aelseif (j==k&&i==j)

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

                    temp=0;                   %thth

0 M7 u$ W# J% B, H/ t. F
for l=1:30

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

; h; C& `! Q4 dend

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

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

( N9 k/ c- [& I2 \* a( V: C' @
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);

) K* a6 P/ C+ F  z1 `  Ielseif (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);

: b3 L+ G3 ], u
end

9 `8 V0 }* t" aend

# \5 s5 B1 f8 g, ]. A) a1 t
end

4 b, V, O: ~3 @* d( H
end
  ~3 f9 W; F+ C9 B! s6 K%ÖÁ´ËÒÑÐγɣ¨13-42£¬13-42£©ºÍ£¨42-72£¬43-72£©

  }4 U3 ^9 G9 ?5 v2 u( x$ L3 vfor i=1:30

& B' t/ k  i5 x: N& S) M
for j=1:30

3 `3 f# v, M( T2 I/ ufor k=1:30

* t2 x0 r  V' N9 n- ?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

5 Y  |  _) ?  F( ]$ ?
elseif (j==k&&i==j)

                    temp=0;                   %thV

* ]: A. z2 |3 b) a8 u( D* K
for l=1:30

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

1 K# ~6 X. L+ I
end

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

. e9 Z2 S; m7 n7 W0 E- \4 e6 t& d
elseif (j==i)

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

! X+ Z2 ~: M+ J8 s
elseif (k==i)

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

2 G! H% K% z, X) ^% P7 r8 Q4 H/ @2 P2 x
end

. T+ Z4 `" }2 q2 _end

2 Q- t! e* R. \& Q, b0 [end

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

' ]0 \& [/ ?! h( u
end. J3 `* e/ z  \5 l; ?
%ÖÁ´ËÒÑÐγɣ¨42-72£¬13-42£©ºÍ£¨13-42£¬43-72£©

; [3 e4 Z2 u4 i* V7 x7 H. a& S
%HhÐγÉÍê±Ï

5 Y$ x8 C: b' _- l+ g. r%Fourth Step: Jg, Hg

    Jg=eye(42,42);

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

    Hg=zeros(72);

) e- X0 j2 M; S- `4 K# J$ E" u%Calculation Jacobian&Hessian matrix END

, t3 S' f. A3 o0 b$ N%%

9 E' i. _3 u. w" n2 D* V$ b%Calculate Newton Iteration Îó²îµü´úÁ¿

/ A! s! Y) }# K9 ]
%Cal LX0-------------------------1

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

( U. Y; K8 t7 Z%Cal LLam-------------------------2

    LLam0=h;

    pferr=max(LLam0);

; f0 M5 o; B4 `! ^
%Cal LMU_MIN-------------------------3

    LMU_MIN0=g-sl-gmin;

2 S8 ?8 d2 l4 f" `4 Z0 q8 ~  |2 V
%Cal LMU_MAX-------------------------4

    LMU_MAX0=g+su-gmax;

1 c0 f9 y, |4 l! H. r3 c2 b" p%Cal Lsl-------------------------5

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

. C& ]" w4 Z) F9 x+ G
%CAl Lsu-------------------------6

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

# n  g: O) [) k. [2 j%Calculate Newton Iteration Îó²îµü´úÁ¿ END!!!

8 t' W: K  ^6 h- h. s%%

$ ~2 w' Z7 a, d) T7 v
%Calculate Newton Iteration ·ÂÉäÐÞÕýÁ¿

- i+ S7 A! ~( J: W( j. q! T
%1st Step: ·ÂÉäÐÞÕýÁ¿delXaf,·ÂÉäÐÞÕýÁ¿delLamaf

2 U  J& e. O7 O2 V4 F; @%S1:H

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

# ]- J' c2 Z3 J* sfor i=1:60

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

' p; i1 J$ o' d( X6 N2 j) o+ Gend

     tempgMUg1=0;

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

9 t& P1 Q* \, Rfor i=1:42

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

4 g$ M( G) @8 p4 y# C
end

     tempgMUg2=0;

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

8 f5 h" A2 g- F9 t/ |% ^4 e- ~  q2 E
for i=1:42

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

' E$ Y. I- O' Y! v9 t* S9 lend

     H=Hf-temp+tempgMUg1+tempgMUg2;

8 h( z* Y# Y! e$ M  _. \
%S2:·ÂÉäKESAaf

     tempgMUg1=0;

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

* g- a( m1 F9 {7 S
for i=1:42

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

0 g6 q& z* p6 Z8 ^; gend

     tempgMUg1=Jg*tempgMUg1;

     tempgMUg2=0;

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

+ C% @" v  P3 B+ l
for i=1:42

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

) r/ |' q" \1 Gend

     tempgMUg2=Jg*tempgMUg2;

     KESAaf=LX0+tempgMUg1-tempgMUg2;

6 D2 z% A- s+ M$ r/ z( n' G%S3:JACOBIAN

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

9 \+ A& q, E: ^9 k: @' [8 S* s
%S4£ºRESULT

     RESULT=-[KESAaf;LLam0];

3 i7 [. N8 b# P* b3 O, d5 ~
%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

* }; j2 W& X' a2 H" ]8 \

' N  S) s# E% l' k' [
%2nd Step: ·ÅÉäËɳڱäÁ¿ÐÞÕýÁ¿delslaf, delsuaf

6 G6 t% o) J& [9 Q) Z8 K  Q%S1: delslaf

     delslaf=Jg'*delXaf+LMU_MIN0;

+ W8 C# |3 g5 v- f; T9 o9 c* W%S2: delsuaf

     delsuaf=-Jg'*delXaf-LMU_MAX0;

& ~! e5 n1 e! v0 G1 e& W/ E' r) g

. z3 x8 M$ s7 ^1 v, A" L
%3rd Step: ·ÅÉäÀ ¸ñÀÊÈÕ³Ë×ÓÐÞÕýÁ¿delMU_MINaf,delMU_MAXaf

' P- n, B$ k7 Z/ A; D%S1: delMU_MINaf

     temp=0;

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

' X3 O/ B5 I$ Dfor i=1:42

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

/ c: y& ]8 {" M/ o; g
end

     temp(13)=0;

     delMU_MINaf=temp;

, L1 P$ l% F7 x; A& ^
%S2: delMU_MAXaf

     temp=0;

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

* {+ r, N3 n8 ^+ ^" l6 pfor i=1:42

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

' D% k% I" I/ p# I, ^# a
end

     temp(13)=0;

     delMU_MAXaf=temp;

/ e" F7 d9 e$ W! B1 y: o; \- q  W" S

, q* G2 f% b& n( o%Calculate Newton Iteration ·ÂÉäÐÞÕýÁ¿ END!!!!!!!!!!!!!!

5 E/ Y6 j: K* g# B

3 U& T2 A' N9 `3 K2 j# d%%

* p8 h7 C! `% ]%¼ÆËã·ÂÉäÔ Ê¼²½³¤ºÍ¶Ôż²½³¤ STEPpaf£¬ STEPdaf

/ @8 O! X6 g) a! M
%S1: STEPpaf

6 j7 {0 }# V8 E
for i=1:42

$ d7 s# K. ]0 @' c6 P; ~( j, |7 _, Q0 l
if (delslaf(i)~=0)

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

0 P" a5 w0 w2 ~! k" W
else

            temp1=Inf;

, x5 ?9 D, f! ]: T, f& [
end

" R2 L( ?% Z, R3 Q) Lif (delsuaf(i)~=0)

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

6 `6 n, P# [( p2 Y* H  ?else

            temp2=Inf;

0 }$ ~7 n& b5 V0 t: y5 z6 V
end

! Z8 @5 e' \* f. ?; Z7 ?if temp1<temp2

            min=temp1;

. J2 L* k! s) [6 v5 n$ K
else

            min=temp2;

- K; A2 a' U% J! {) l2 W7 `; Gend

2 g" y: t* z- x; uif min>1

            min=1;

! i; H2 d) H) V0 a, {' T8 O7 Vend

% f! c0 s. b9 M; l5 w! I" D! s- P
end

    STEPpaf=0.9995*min;

- `: K6 g% {6 n
%S2: STEPdaf

$ R& r: ^9 a  b# u. D# R3 z, w8 G9 o# _for i=1:42

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

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

3 `. i" \1 X/ E' b6 t! l
if temp1<temp2

            min=temp1;

: I' I3 C  v7 y$ S$ S( Telse

            min=temp2;

' |( G, n$ R; W; Y, I3 z: ~
end

+ V( S( y) d4 Wif (i==13)

            min=0;

  t" A/ o0 w9 ~; cend

; c5 n: K8 Y. D5 W& D0 N* J" a9 T6 x! hif min>1

            min=1;

9 J( `. Y0 X( n" W7 N$ Q  d3 p8 g) s+ N
end

* r5 H$ e% ^6 D" a3 x5 y* }# N" eend

    STEPdaf=0.9995*min;

4 z' M' o$ `; }8 R%&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!!!!!

/ t# A  l3 a; \% r

3 ?4 E1 m, f' b4 _
%&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);

$ F) G7 ^0 t+ h! h
if (ROUaf/ROU)^2<0.2

       SIGMA=(ROUaf/ROU)^2;

3 T5 o, C+ f( c% O4 O1 g' zelse

       SIGMA=0.2;

* p1 [. ?# @+ m# i
end

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

* q: J% B  G1 P5 o: U+ 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;

" g, c8 M& t$ H, c5 _* U" }! ?( R
%Calculate Newton Iteration &ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;

6 {: ^$ H2 g/ \" d; c: [
%1st Step: &ETH;&pound;&Otilde;&yacute;&ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;delX,&ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;delLam

9 ^0 D6 K. G- R# e6 p. V
%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;

" Y- r  ]8 u+ c2 Jfor i=1:60

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

0 N& u' Z6 [/ `, e1 _, R$ D- s" m
end

     tempgMUg1=0;

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

. V2 _% N( K8 ~# \2 V! }9 B
for i=1:42

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

* d$ L$ O  p1 l% h# |- Fend

     tempgMUg2=0;

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

0 }7 M  [! z9 w" _- i( wfor i=1:42

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

) |* |# e0 i  P7 s, d. N' Y
end

     H=Hf-temp+tempgMUg1+tempgMUg2;

% C- v+ m9 x4 U  c%S2:&ETH;&pound;&Otilde;&yacute;KESA

     tempgMUg1=0;

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

. Y+ E5 ~, C1 e) Y6 S( A/ Pfor i=1:42

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

) V" [* ~  ^' g% o8 v3 p
end

     tempgMUg1=Jg*tempgMUg1;

     tempgMUg2=0;

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

( y, r, m7 N. d0 j- j
for i=1:42

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

/ w; }4 |% ?( ?2 L4 d$ w
end

     tempgMUg2=Jg*tempgMUg2;

     KESA=LX0+tempgMUg1-tempgMUg2;

; p9 X" f- u% W, P' c%S3:JACOBIAN

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

: ?1 |/ R% R1 |; w- t%S4&pound;&ordm;RESULT

     RESULT=-[KESA;LLam0];

8 P( ]6 c8 A0 ~) w! e%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

* ~$ P+ h, u9 k. h9 K9 u( d  _, C

7 T5 H, Z* {; T8 G8 }, \%2nd Step: &ETH;&pound;&Otilde;&yacute;&Euml;&Eacute;&sup3;&Uacute;±&auml;&Aacute;&iquest;&ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;delsl, delsu

: i9 o8 i- [  J& z7 l4 U  ~9 @
%S1: delsl

     delsl=Jg'*delX+LMU_MIN0;

) n- ]7 P2 }$ L6 Q7 z5 H1 [
%S2: delsu

     delsu=-Jg'*delX-LMU_MAX0;

2 Z! n% Y) p- A3 l

( _: S& q5 U9 l. I2 J  Y%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

, f$ p+ R( w+ K; N%S1: delMU_MIN

     temp=0;

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

- d2 @1 f  p2 c" k
for i=1:42

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

6 B( z6 X" \  K  G$ b, U" Jend

     temp(13)=0;

     delMU_MIN=temp;

7 L; z" c# F) Z$ j& C5 P  a
%S2: delMU_MAX

     temp=0;

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

3 B; `( d1 k1 f4 I
for i=1:42

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

+ T' N% X3 c3 E6 u
end

     temp(13)=0;

     delMU_MAX=temp;

+ N6 ^4 @9 L8 W1 d+ I! Y

! U( N' a: i& p1 D4 B
%Calculate Newton Iteration &ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest; END!!!!!!!!!!!!!!

- ?4 n2 y1 C' w3 S
%%

! r% s- u2 F( D$ n& E%&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

2 U$ G- E& A, K8 i  {2 G6 M! w* e%S1: STEPp

2 J0 M% w" D2 m3 zfor i=1:42

+ z* c* h6 V. Y* H* @) Q- h* w
if (delslaf(i)~=0)

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

% l9 Y; ?8 f4 gelse

            temp1=Inf;

; i+ F. b! k8 x& i3 v& i
end

& I3 R7 i! A/ U+ Q( `1 {, M- mif (delsuaf(i)~=0)

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

) ?9 T- @7 |9 yelse

            temp2=Inf;

7 w1 X2 z# y  D& P5 h, Hend

  A- `! O9 c' R/ \* a7 mif temp1<temp2

            min=temp1;

( A4 ]9 b6 B" `$ kelse

            min=temp2;

8 J* g' l: l+ h. d7 lend

( d7 H2 B& ]0 Qif min>1

            min=1;

3 A" X" M, q. v' a
end

# F& J3 S8 I( H; Z# `/ Y+ P' V
end

    STEPp=0.9995*min;

; E" q3 W5 W9 E% v5 S
%S2: STEPd

* y( T6 ~3 N, }! |! Vfor i=1:42

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

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

5 E  ~+ m  \1 c
if temp1<temp2

            min=temp1;

. Z( e* Y& K6 ]9 Kelse

            min=temp2;

& @; U1 X' c# G* X; E. }% O0 [& T7 tend

3 X. d9 U5 V% \: N! S2 s6 Zif (i==13)

            min=0;

) E* Y7 G) y- u0 H
end

" Z3 f3 j( ?/ O$ w8 U" Jif min>1

            min=1;

5 j$ j$ N, A( F. r, M, Cend

! M' i. o( W; ^' R  X( q  ]  v' Pend

    STEPd=0.9995*min;

2 Q8 b/ q" {# `/ L1 U  M
%&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);

; b8 z4 U: @4 @8 w
%&cedil;ü&ETH;&Acirc;&Ocirc; &Ecirc;&frac14;±&auml;&Aacute;&iquest;&ordm;&Iacute;&para;&Ocirc;&Aring;&frac14;±&auml;&Aacute;&iquest; END!!!

/ i0 z. M# F3 k# c6 U9 H

& _0 r  Y+ `- g5 }. ]% u%&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));

: t2 z2 W" _# P; U
%%

8 m9 K: A- x( E" ^2 }4 [8 k; D
%&Aring;&ETH;&para;&Iuml;&Ecirc;&Ccedil;·&ntilde;&sup3;&not;&sup1;&yacute;&micro;ü&acute;ú&acute;&Icirc;&Ecirc;&yacute;

    ik=ik+1;

! L# _% r5 M: {& K) C
if ik>IterNumMAX

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

/ ]8 P* F/ J5 P7 o
break;

8 N5 h6 ~& e4 v
end

% L0 V$ _, }  _1 p1 F" u

end

et = etime(clock, t1);

%%

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

/ l7 G( Z7 b3 k+ w

if (ik<=IterNumMAX)

% _1 U# m% Y" {- }  \
%%

# A$ \$ T# u0 S. ~( y
%&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

$ Z8 d, X- r* \

max=-10;min=10;

for i=1:ik

    temp=LLam0(i);

2 B, M8 \6 b8 p( p) Q3 {- M
if (temp>max) max=temp;

8 |  V. Z5 m, V5 ~, vend

  {+ X4 k. P; X. H1 s  ?/ G& ]
if (temp<min) min=temp;

7 I% k. J1 Z7 f  b% L) D, M( M
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('  -------  -------');

2 i) C) A6 m9 i3 R! k9 T
for i = 1:30

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

5 d8 k3 M" Z# l" P! h6 B) ]7 m3 Fif (i<=6)

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

0 Z4 h# S+ J; i, z7 t
else

            fprintf('       -         -  ');

5 J5 w0 K$ H5 j( ^) o6 I) yend

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

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

, T" l2 n4 H+ b) y/ {
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');

* @8 u" |3 z. M7 v4 ^7 Z: L: M
end

- c. \- Q" v9 ?+ Q3 g( Oend

' T8 a# `* q4 s9 z) E
%%

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