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

! ^: r1 u$ h3 I%%

5 r; |% j. Z. a$ V  P. B" K
%Calcute h,g matrix

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

        errArr=[errArr;ROU;];

    SIGMA=0;

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

0 l" U" z, f' E8 G9 S
for i=1:30

        temp=0;

% Y- |# Z! t$ [6 a  K* \. a/ qfor j=1:30

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

/ j+ r* L  x8 G: T7 t/ A( Q
end

, D+ D3 M) K4 D% ^: ]" t! Q/ ?/ Sif (i>6)

            tPg=0;

6 u$ e0 H8 @3 Z% b8 Eelse

            tPg=Pg(i);

- {( ]6 a& w2 a8 `, K/ [1 |
end

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

" q3 s" a2 c1 L  h& F! B( k  E0 [2 [
end

, o# y. I& K" ~5 ~4 y& k4 p! z0 }; s  gfor i=1:30

        temp=0;

2 p% G2 N& t0 Z* [0 `for j=1:30

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

1 i' a6 \: Y' Z2 n4 ^) o: Qend

6 e4 x0 y8 i' T0 }9 n6 J
if (i>6)

            tQg=0;

, s% `0 S5 |; k1 \8 _: O% E2 z
else

            tQg=Qg(i);

& B7 f0 M9 G0 Z8 T
end

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

" o3 ?% G, _* |+ G
end5 b. R4 t: i' u0 u* n8 e7 U; b5 k7 C
% Cal  h END

2 \# F$ `! t9 z' y" @

2 V7 F1 D" C. {- l% f+ D
for i=1:6

        g(i)=Pg(i);

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

. ~$ H- {. j- _- m" t& Tend

8 y0 o8 }7 d  a+ ^' K4 Z2 M. a
for i=1:30

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

1 M% M. _6 X$ D( ]8 d1 @
end
9 L6 k. b' C1 k, r" U% N5 j% Cal  g  END

/ t: t% c  G  X5 u%Calcute h,g matrix END

: x$ ?6 y/ q5 g7 N%%

. o1 u7 [4 Q3 Q6 y/ p, v' D
%Calculate Jacobian&Hessian matix

5 E" N) ^% ]# v! x) |%First Step: Jf,Hf

  T0 M/ U# ?4 o$ K0 A1 k9 r
for i=1:6

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

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

9 {1 w1 Z9 V0 Y2 x! Z6 d
end

/ s' M7 i2 x, o9 l3 L$ r%Second Step: Jh, hΪµÈʽԼÊø

2 R4 ~, W0 |. W3 G# Y# @7 Pfor i=1:6         %Ç°6ÐжÔPgÇóµ¼£¬ÓÉ´ËÒÑÇó³ö

        Jh(i,i)=1;

$ `& ~7 G8 K* H5 _7 @: ^
end

# @% L- @  b+ p* m; w$ q: ufor i=7:12        %7-12ÐжÔQgÇóµ¼£¬ÓÉ´ËÒÑÇó³ö

        Jh(i,i+24)=1;

# O3 F" Z1 h6 }4 `( ]8 k
end

/ n$ k8 `3 `8 C1 M" u
for i=1:30       %ÐγÉ13-42ÐеÄ1-60ÁÐ

- K: B! z# R+ j! J) }
for j=1:30

            tempVp=0;

            tempVq=0;

# j0 x+ u" d4 i( f! e' Qif (j==i)

9 w  D: Y; g( I# a9 V5 y, {# B
for 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));

' [: h0 t3 }/ h8 m3 U- W
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));

* _# Q3 d! g( H9 {& W) d. eelse

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

' c; r' p7 Z! U! \( A3 W
end

% v4 [+ R0 Y5 I3 Nend

2 K! T* {  r. X2 y4 _3 b! N7 f! ?end

1 S" N. G7 k+ ^" O. }
for i=1:30       %ÐγÉ43-72ÐеÄ1-60ÁÐ

+ \, o( y: x& V4 |for j=1:30

            tempVp=0;

            tempVq=0;

' _3 n" k- Q0 f, E4 g! v; jif (j==i)

0 p/ r' B% i# j/ I  U; G1 Y/ w3 W" Vfor 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));

; v4 S7 A+ V- C
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;

& m% f  \! ?9 |( Gelse

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

7 X6 [; |; C! Z) P+ g/ Send

: B3 q7 R9 k* I/ V+ a  }+ \
end

1 b8 ?/ m: z5 B8 g! X3 R) g
end

6 g9 ~& O% h! A+ N" {
%Third Step: Hh

. _# Q# @! g! _$ m5 _%Óй¦²¿·Ö

3 m1 P$ h1 ~3 Y( @, a
for i=1:30

* s4 p; R" e9 g6 ?* Gfor j=1:30

$ l& @  L0 `# s$ b
for k=j:30

. c: A! g9 G, W* O7 [& Z& Z. [
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

6 h, V' ]8 ~( R$ E* X1 S7 G4 Uelseif (j==k&&i==j)

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

                    temp=0;                   %thth

$ g9 Y2 Q% |' z8 a# d3 S' Vfor l=1:30

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

2 U& N9 Z- k& Cend

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

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

* Y% d& I. x0 ~% G
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);

: U' Q/ v  g0 G7 K' P
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);

/ w2 L% r8 c0 iend

/ w6 R) b  Y# d  V( \& mend

) k6 o4 j  i* x, X+ V
end

, y1 }8 M: ?* m7 C9 y: u$ l1 Vend; W3 ?4 g7 H& Q  C0 W1 O2 r# c( O
%ÖÁ´ËÒÑÐγɣ¨13-42£¬13-42£©ºÍ£¨42-72£¬43-72£©

" T- y2 y8 L- y4 Zfor i=1:30

8 z7 |: r- N: s+ l8 i- I& d3 ^
for j=1:30

! L, N! f+ Q8 d5 m8 ]0 P
for k=1:30

8 c) u& E& u2 _5 {; }
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

0 b5 x# c9 r- {
elseif (j==k&&i==j)

                    temp=0;                   %thV

: w6 |" l7 q' }# d6 L7 U5 I" u$ x) zfor l=1:30

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

3 q% y. f# I: {/ i- D1 Jend

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

+ M  L# H1 I. q
elseif (j==i)

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

0 J- L) i6 N$ l9 z' k) s! felseif (k==i)

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

  E) N# U# }3 _0 F4 K/ h; }8 \
end

/ R' O' C: }* N& a# A
end

* q7 _+ v: R5 p$ I
end

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

3 m: ?3 V% p5 U5 J2 {end$ q! ]9 r* S" d9 x! z2 i
%ÖÁ´ËÒÑÐγɣ¨42-72£¬13-42£©ºÍ£¨13-42£¬43-72£©

, P) r+ [+ k0 r! [* F%ÎÞ¹¦²¿·Ö

, @1 b9 j( M  `" Qfor i=1:30

5 ^* g0 Y# d: Z+ ~
for j=1:30

! F$ a9 G  f* S2 Pfor k=j:30

  v7 h" I: m1 Q) ?/ [  b& Z3 t' Sif (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

( w) P/ r! z( B9 }3 Belseif (j==k&&i==j)

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

                    temp=0;                   %thth

, F% s7 E: r. s4 ?) I% ifor l=1:30

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

- K$ d& d# {( K: A: e
end

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

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

, C+ i( C  u  e' K9 X- A2 kelseif (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);

& U7 q, r. f0 _: C
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);

: K+ U* G3 G9 \4 c
end

% d3 q) l2 ]; W4 j) K5 z5 A
end

# N0 L! ]* F1 c: Q' Cend

" m" v( Q+ y1 Z. }, m7 @end
- g5 }4 h( u/ p) e%ÖÁ´ËÒÑÐγɣ¨13-42£¬13-42£©ºÍ£¨42-72£¬43-72£©

' T! P4 v% C, r7 T2 F2 R
for i=1:30

. V) ^0 B+ s9 \: [
for j=1:30

6 o% X% a( |0 e5 n! A( Z
for k=1:30

" p) |: Y2 [8 E$ {" f- t
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

$ {. T6 E0 O. R/ S3 }; m0 c
elseif (j==k&&i==j)

                    temp=0;                   %thV

5 P9 @  Z1 B. k5 ?4 H
for l=1:30

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

* h- [/ O1 R  Q& A- U, S: o9 i2 Gend

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

+ }8 t1 ]  m+ ~: F% k- B3 e( {
elseif (j==i)

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

  `3 b, g  h3 A9 ?; W! relseif (k==i)

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

! e0 V4 O( I  j# F; ]
end

& @6 @# E0 A$ ^- \' _end

0 o$ H* y0 W2 v. fend

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

2 K' _% ]! Z9 |. d4 z
end
4 q9 @0 L: m" l" Q( i%ÖÁ´ËÒÑÐγɣ¨42-72£¬13-42£©ºÍ£¨13-42£¬43-72£©

& v! b" p$ K' F/ K( H%HhÐγÉÍê±Ï

- K. ]( L! f. K* ?3 I$ ?
%Fourth Step: Jg, Hg

    Jg=eye(42,42);

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

    Hg=zeros(72);

5 S2 Y: ?: s5 s2 ^) }) @%Calculation Jacobian&Hessian matrix END

  i) G3 W. I" U1 m1 I4 m
%%

& y8 d+ s- d% T  w1 o. Y- F2 G. w( u
%Calculate Newton Iteration Îó²îµü´úÁ¿

- C, ^: R/ R  q  P! o) x
%Cal LX0-------------------------1

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

- A" `9 g+ G4 I1 w& O3 V* F
%Cal LLam-------------------------2

    LLam0=h;

    pferr=max(LLam0);

! {3 d, g6 w0 L' g7 j' i: _
%Cal LMU_MIN-------------------------3

    LMU_MIN0=g-sl-gmin;

+ Z' \. {  ?6 m/ ~. s) [) x
%Cal LMU_MAX-------------------------4

    LMU_MAX0=g+su-gmax;

) i: O( B7 ]4 I%Cal Lsl-------------------------5

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

- `" O: [6 g& _. m6 z, U
%CAl Lsu-------------------------6

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

; J: E. j8 ~( b7 ]; R% s
%Calculate Newton Iteration Îó²îµü´úÁ¿ END!!!

0 x; I( y8 K" S' X" s0 g8 I
%%

" V; h1 w$ V% w! x1 l2 a
%Calculate Newton Iteration ·ÂÉäÐÞÕýÁ¿

* x0 m% C# k; b9 r%1st Step: ·ÂÉäÐÞÕýÁ¿delXaf,·ÂÉäÐÞÕýÁ¿delLamaf

1 d) Y0 d$ `  W$ [; Y" \, Q" [%S1:H

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

" M  h" ^# H: u. [: Ufor i=1:60

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

5 k% O7 }5 r+ Y8 S- W( [4 |end

     tempgMUg1=0;

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

2 i0 X& M0 N8 Q$ Mfor i=1:42

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

& b" }/ w4 f" v4 h& K/ X. H0 r, Y
end

     tempgMUg2=0;

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

; e* W1 f1 ?, b, Yfor i=1:42

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

4 S( B: P3 [: X* ]9 E' r6 B0 Dend

     H=Hf-temp+tempgMUg1+tempgMUg2;

5 I, n2 E6 R* o% Y- Y
%S2:·ÂÉäKESAaf

     tempgMUg1=0;

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

9 \+ J. L5 k8 L
for i=1:42

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

5 Y1 x( X2 F3 c/ U5 k2 l0 Y' {4 Eend

     tempgMUg1=Jg*tempgMUg1;

     tempgMUg2=0;

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

* f2 t. l! l! i  T' Wfor i=1:42

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

; E$ q7 _; _6 s5 d/ U- J3 Z5 Nend

     tempgMUg2=Jg*tempgMUg2;

     KESAaf=LX0+tempgMUg1-tempgMUg2;

; A5 U, _6 h2 A9 w  C- a
%S3:JACOBIAN

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

$ t5 t/ a$ K% Z$ f* `$ D
%S4£ºRESULT

     RESULT=-[KESAaf;LLam0];

' O: i. K& T! D( Y: Z" X5 ]* O%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

3 z+ {# m6 G. Z, G' s5 g5 B

' F# ~- p; V9 {# f: ^5 x# h* q- l
%2nd Step: ·ÅÉäËɳڱäÁ¿ÐÞÕýÁ¿delslaf, delsuaf

, y# J9 ~% R5 H%S1: delslaf

     delslaf=Jg'*delXaf+LMU_MIN0;

5 ^  d. O3 I% F9 j. u& f9 k%S2: delsuaf

     delsuaf=-Jg'*delXaf-LMU_MAX0;

* m6 o+ J. }3 z7 E6 Q2 Q; s

- P5 I. G# B( V. W& c%3rd Step: ·ÅÉäÀ ¸ñÀÊÈÕ³Ë×ÓÐÞÕýÁ¿delMU_MINaf,delMU_MAXaf

' [! B' t2 u2 L7 J) l4 ]3 e; j%S1: delMU_MINaf

     temp=0;

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

/ d) j% A5 u; F# V5 i# l
for i=1:42

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

4 [- M$ i. k$ Y1 z! eend

     temp(13)=0;

     delMU_MINaf=temp;

, V$ ^/ ?! v6 g3 p1 X2 h%S2: delMU_MAXaf

     temp=0;

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

9 Q$ Q  s; }- l/ Z) G  Z8 tfor i=1:42

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

; L, R% [8 S! B# F, d9 kend

     temp(13)=0;

     delMU_MAXaf=temp;

) i2 p4 G' i* `& k' o2 j( z- T

! |0 U2 E1 |4 l! B% @%Calculate Newton Iteration ·ÂÉäÐÞÕýÁ¿ END!!!!!!!!!!!!!!

+ U$ ~. S  y* H) S4 o9 A( R! @0 L

* x8 i9 H  o3 y  F2 Z%%

- d5 q2 L! h  ^8 x- I. E7 K4 W
%¼ÆËã·ÂÉäÔ Ê¼²½³¤ºÍ¶Ôż²½³¤ STEPpaf£¬ STEPdaf

( d+ ~3 S4 k! H$ o%S1: STEPpaf

. ~% F% [! |2 t+ vfor i=1:42

0 ?; Z  \7 R# f4 kif (delslaf(i)~=0)

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

/ H6 G; c% v; melse

            temp1=Inf;

: f: k! x0 G, I2 f) V. d; z0 dend

6 V/ a( K3 Y: d: ]. I, e! gif (delsuaf(i)~=0)

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

/ L, s) h8 m! N: Z9 qelse

            temp2=Inf;

, ?( a- f5 e: t2 F7 N& ]+ g. l
end

* }2 C' I4 y2 R% e( h% v7 cif temp1<temp2

            min=temp1;

+ g* t$ x$ Z8 e* I* r
else

            min=temp2;

, G6 s) H9 Y! m. t( _end

1 V& @4 q0 ?" `, {& ]9 j
if min>1

            min=1;

1 o- }4 Y. T, i* R# [. ?; Tend

$ O! c7 J9 Y& o$ |end

    STEPpaf=0.9995*min;

# I, ~4 s! O% X5 b9 D1 t1 N0 l% K
%S2: STEPdaf

8 E, e  a+ Z3 dfor i=1:42

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

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

" V* w* U) w4 [1 E4 Nif temp1<temp2

            min=temp1;

# J, n, H7 S* L1 C* e& j
else

            min=temp2;

# R4 ~# X' j& \end

; d8 p+ N9 c  R* O( D
if (i==13)

            min=0;

) M( _: i* O) k! [0 R1 Cend

5 Q* ~1 s5 M5 X4 B, S7 Bif min>1

            min=1;

; f% N: I3 B2 G- P
end

4 I; `, _/ i+ Xend

    STEPdaf=0.9995*min;

3 r4 z/ B1 o* @+ E& ^* q
%&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!!!!!

8 O: h* c: S" p! H' o8 K

# [( {7 F4 v8 I1 i- O% A6 r%&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);

0 W& W: ~; [+ X% ^
if (ROUaf/ROU)^2<0.2

       SIGMA=(ROUaf/ROU)^2;

& }9 ~1 @% Y3 L, C- u* F, u- g
else

       SIGMA=0.2;

  o. q8 \7 H  V0 a: ^
end

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

& Q; Y; a$ E* x7 K1 f%&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;

3 t2 }2 y, o/ E
%Calculate Newton Iteration &ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;

, Q: `/ B7 {4 Y%1st Step: &ETH;&pound;&Otilde;&yacute;&ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;delX,&ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;delLam

- `6 W& N+ ~% Y5 _" L0 X%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;

1 o) y# S" G  c/ Nfor i=1:60

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

& A% l: O) }! Z" oend

     tempgMUg1=0;

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

/ `7 Z: u+ q' Q. tfor i=1:42

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

6 K0 {: @& K$ R  n) e8 ~% y
end

     tempgMUg2=0;

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

0 N4 ?' p" [! j8 C" e( [
for i=1:42

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

! F6 ]2 u" K' N' C$ A; zend

     H=Hf-temp+tempgMUg1+tempgMUg2;

1 Y3 p# n; `- g! b* }
%S2:&ETH;&pound;&Otilde;&yacute;KESA

     tempgMUg1=0;

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

6 i4 r3 }2 u. H. A; r' U' m' p8 a2 |
for i=1:42

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

% _. ?2 n; C& Q+ b/ |2 n$ L  p7 y# qend

     tempgMUg1=Jg*tempgMUg1;

     tempgMUg2=0;

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

% _0 R+ U! h, R* t/ p  ?
for i=1:42

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

8 D: t% q) n7 qend

     tempgMUg2=Jg*tempgMUg2;

     KESA=LX0+tempgMUg1-tempgMUg2;

& s! L0 d, A% b) l6 ]- u5 ]
%S3:JACOBIAN

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

7 n, t* |9 Y7 d3 W5 ^- ?5 O4 t! }%S4&pound;&ordm;RESULT

     RESULT=-[KESA;LLam0];

. Q' u! x2 t" N2 R1 W* y1 z
%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

! m: E6 x3 F( x

7 K8 W; z8 I$ l* B%2nd Step: &ETH;&pound;&Otilde;&yacute;&Euml;&Eacute;&sup3;&Uacute;±&auml;&Aacute;&iquest;&ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;delsl, delsu

; c: ]8 Y2 ]( `! D! l%S1: delsl

     delsl=Jg'*delX+LMU_MIN0;

$ P! \' [# @' L/ M% d7 l%S2: delsu

     delsu=-Jg'*delX-LMU_MAX0;

% O! a4 f6 w1 N7 T$ @

8 i( @; _4 M( }2 y' y( ], b3 T%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

- t$ u, V2 @" F0 c  s%S1: delMU_MIN

     temp=0;

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

* S7 N1 T! F6 Y& k6 A: Dfor i=1:42

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

3 R9 S$ K: _( k
end

     temp(13)=0;

     delMU_MIN=temp;

/ a- X$ l" R. ~% V
%S2: delMU_MAX

     temp=0;

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

3 i/ i! Q1 ^( J) `2 Z
for i=1:42

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

* u9 U; w, A) j0 _" A/ b  ^
end

     temp(13)=0;

     delMU_MAX=temp;

6 C% p6 e3 G* v! ~: r" M  h

7 q7 B$ S2 H4 @
%Calculate Newton Iteration &ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest; END!!!!!!!!!!!!!!

+ i7 |5 o( o- Y' A+ a%%

% ?) g2 H+ ?  I( T%&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

8 R: [  v* C" {' H) ?5 |+ B( N
%S1: STEPp

1 h" M8 q! ?) n, F" @for i=1:42

& d7 n0 {/ b* A& r: k9 M, jif (delslaf(i)~=0)

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

( k- y5 r- L. t; l! p) b
else

            temp1=Inf;

" X, ?% x, ?8 S! x& M+ J7 Q8 nend

6 g$ p, q+ r( A2 I" N! A4 q
if (delsuaf(i)~=0)

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

) s; F2 l8 I$ \7 @else

            temp2=Inf;

( V) H) N$ J; cend

. E6 M2 l0 ^- {2 M8 |* T# z; V# {if temp1<temp2

            min=temp1;

& I. w5 j8 V1 y' Y( ?else

            min=temp2;

. q6 Q$ I4 {9 c) O: N
end

1 |! s9 F  @: e( T9 U# b6 Vif min>1

            min=1;

: l/ o  ?* o' g/ t0 o* }+ C* Vend

2 V. _4 ~) z% ?. }  e3 e+ rend

    STEPp=0.9995*min;

4 u5 l  P4 L9 |
%S2: STEPd

0 w! i+ ~1 m9 r1 n) u& p  c+ Tfor i=1:42

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

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

% N% h$ U' P( W# K( kif temp1<temp2

            min=temp1;

3 y& `8 P5 N8 H6 f( Q( K) V
else

            min=temp2;

1 e9 L' s5 f! F# Y' Fend

0 i9 q# q+ o, b4 o: T" x9 zif (i==13)

            min=0;

, ?) z1 B8 n/ v+ D
end

) q) z  R. ^$ h/ v
if min>1

            min=1;

" |: T: r' {6 F# N
end

' `2 y5 U. a$ s
end

    STEPd=0.9995*min;

3 F+ A: j# I+ 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);

' h' t$ r- t* c7 m4 w% m# Q4 m%&cedil;ü&ETH;&Acirc;&Ocirc; &Ecirc;&frac14;±&auml;&Aacute;&iquest;&ordm;&Iacute;&para;&Ocirc;&Aring;&frac14;±&auml;&Aacute;&iquest; END!!!

& N/ K6 W& K! y+ \/ A

. N3 \; q1 f, g9 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));

8 ]; m! C# n3 w- [" O7 R$ U4 \%%

& R2 {) a5 F9 D* X, J  a. i4 e
%&Aring;&ETH;&para;&Iuml;&Ecirc;&Ccedil;·&ntilde;&sup3;&not;&sup1;&yacute;&micro;ü&acute;ú&acute;&Icirc;&Ecirc;&yacute;

    ik=ik+1;

+ M2 n4 W3 O% L5 y+ d3 Hif ik>IterNumMAX

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

9 X* U, G/ P8 w2 o' W$ B- {  d
break;

* |; D+ O: g$ ~1 p5 k
end

( ]! [" C2 P& {2 X4 O

end

et = etime(clock, t1);

%%

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

$ o3 x; W1 U! G; S# H5 J0 l

if (ik<=IterNumMAX)

7 }5 ~/ M1 Q9 z  g- y$ I
%%

9 W" V5 P/ |7 a. f# `5 b! F" 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

0 P: m& a: F. I  E& Z% B

max=-10;min=10;

for i=1:ik

    temp=LLam0(i);

3 g" Y* @* k% k! I2 Gif (temp>max) max=temp;

. w# l5 N; I3 I; O* S
end

( U8 R0 N0 J4 D; ^$ O
if (temp<min) min=temp;

6 \' g$ B& S3 Q* S( Q1 q, Qend

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

: U2 w7 F9 D0 K+ `& c
for i = 1:30

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

# j/ z" z. f/ f3 E: @- U8 b6 i
if (i<=6)

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

1 v! G5 t$ {+ x" p' ?+ ^4 a
else

            fprintf('       -         -  ');

! l% H1 K% d& [4 b! dend

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

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

  _: c5 H. R& L/ b. R& J3 `, r8 ~2 \
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');

4 V$ f; Q) _. B- h0 J
end

. r6 q- \) T: `6 N' f5 qend

: F% P! k5 Q; n: A%%

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