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

% C3 B) V$ n% ]" S, p$ `
%%

, ?' q; a. l8 p. w9 a7 V) l%Calcute h,g matrix

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

        errArr=[errArr;ROU;];

    SIGMA=0;

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

6 ~6 Z4 E9 h5 ]& j$ v" T
for i=1:30

        temp=0;

! M! N& |9 @/ m3 D9 n
for j=1:30

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

% G( h  v+ ^8 L
end

6 Y. ?# K: x- H& c# S+ G9 Q3 |
if (i>6)

            tPg=0;

. E- V( M+ I) y1 {9 i6 Xelse

            tPg=Pg(i);

- ^! f/ F+ v, ?# j5 I7 R: a0 `end

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

( k4 ~# t, C: V) c8 ?end

6 \0 @6 V  u( I# G/ Ofor i=1:30

        temp=0;

4 n2 O$ d4 X9 s9 S$ a
for j=1:30

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

5 M; ~1 P/ F: [0 i2 b0 ^  @
end

) c' L8 i6 `! {8 _: C6 |if (i>6)

            tQg=0;

( q# ?  L; ~3 C5 d6 W3 A6 eelse

            tQg=Qg(i);

6 D8 S% P; P5 j- ]# }end

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

0 O; L1 [, M1 a( K  {5 F
end. ?( C3 i$ R: M5 R% }- {
% Cal  h END

1 a0 c+ e+ x( s% ^: f( v$ F. ]2 w

; |* f% ~$ e& ]+ g$ W# u- U2 l! h
for i=1:6

        g(i)=Pg(i);

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

! M. n* Y* {) p4 Y
end

, E) x( {+ K7 vfor i=1:30

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

1 |: [) g" m) K2 ?( x% o& ]! q
end
/ y) l8 `* d8 m' K% Cal  g  END

9 t( g9 d/ s, h- p% D  T
%Calcute h,g matrix END

4 p; _# K. D" T0 m%%

" R0 k$ F) `8 M/ U5 t. y%Calculate Jacobian&Hessian matix

5 C% @: _/ J1 N* C%First Step: Jf,Hf

# H0 V) z+ w; n8 b' }7 O& W1 ~
for i=1:6

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

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

4 N3 b+ U0 N: d, B# [6 X  b
end

6 j0 t% p: E5 {2 Q
%Second Step: Jh, hΪµÈʽԼÊø

/ d/ B& b8 M7 Z5 ]% sfor i=1:6         %Ç°6ÐжÔPgÇóµ¼£¬ÓÉ´ËÒÑÇó³ö

        Jh(i,i)=1;

  q: Q. ?# D# p! ?7 x' bend

7 u. Q, p2 i: c! K* V7 e8 D) Y# Nfor i=7:12        %7-12ÐжÔQgÇóµ¼£¬ÓÉ´ËÒÑÇó³ö

        Jh(i,i+24)=1;

( y( l6 j6 C/ _! \, A+ Xend

+ r& H! z- I7 q
for i=1:30       %ÐγÉ13-42ÐеÄ1-60ÁÐ

# ], u& G- {2 i( |2 u& m
for j=1:30

            tempVp=0;

            tempVq=0;

6 h- l( x& ^9 N8 D3 i. Yif (j==i)

% K! v" ]+ p( A, ]$ Q5 }: t
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));

7 V# z+ [4 X% R9 k% Kend

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

, K5 j# P# }5 t4 ]$ q* Jelse

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

. l, ?6 u3 E1 E- I% ~8 Hend

7 [; F* X. `# L4 H: G7 A; K/ p
end

8 b9 W' d) w3 [. M& i! L' \end

* C3 |. w* l$ P
for i=1:30       %ÐγÉ43-72ÐеÄ1-60ÁÐ

9 L: F, h& J/ V! f: p5 m. v, p
for j=1:30

            tempVp=0;

            tempVq=0;

( q% u$ ?$ S. K; p
if (j==i)

5 K' j% i3 J8 F# h7 Hfor 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));

7 A* o: }1 ?' v* D9 X
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;

; f4 f7 [7 O3 {6 ^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));

& g- o+ ~7 s  _4 D4 e
end

& ]$ x8 @. V2 j0 _. rend

* R' `; o  c3 R% z, N3 ]
end

8 s- I# R+ ~1 |7 N+ X8 q$ [1 m
%Third Step: Hh

) m  n! V8 g4 f0 F  i+ p: p0 b
%Óй¦²¿·Ö

+ ^% u6 ]7 c2 y. g1 J3 d0 V
for i=1:30

& ~' f' M/ [) Vfor j=1:30

1 v9 {! {7 g1 K5 j2 d( P$ ^7 N1 C
for k=j:30

! i' _3 Q  u1 B; K0 x
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

( ^4 j: n3 p9 |: e; q: w% A
elseif (j==k&&i==j)

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

                    temp=0;                   %thth

4 t: K' {* Z( V! ifor l=1:30

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

. V9 k6 Q8 `7 L, @& I! mend

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

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

! o- o1 }8 u* r& ?8 ~3 O' W$ p
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);

! l9 ^! d, T$ b$ w
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);

. S7 Z& v6 |: I% D8 s
end

# V( B, p, ~2 \8 q) }
end

4 W" I) W$ r* z  _+ _end

* g4 v1 N! q1 c! F/ q6 mend
( h* s- h8 N6 S. r& |& h! w, F%ÖÁ´ËÒÑÐγɣ¨13-42£¬13-42£©ºÍ£¨42-72£¬43-72£©

+ b2 b1 s. t! _0 v/ K
for i=1:30

  c& g1 t. _3 z& Q
for j=1:30

6 u9 Z9 u* _: E) u. w
for k=1:30

5 b  I6 R5 |/ B5 w, f& Uif (j==k&&i~=j)

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

$ C- f+ B* m, e) v. v
elseif (j==k&&i==j)

                    temp=0;                   %thV

  @: |+ ~! N$ xfor l=1:30

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

( K5 i: |! b2 d, ~3 [end

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

# K, e0 t" ]4 j
elseif (j==i)

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

5 H9 y: ~  }% K# y  }- {' k$ oelseif (k==i)

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

% _7 p" T: O/ X7 _- lend

$ p; j6 ?/ i8 h1 h( \( V# n. h  Yend

( g7 v! O( F! X. v8 N! ~: \end

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

3 `* |, M: @- V  E7 z& i' A- l4 ]end
5 W3 H' \/ m, p5 P' J  m%ÖÁ´ËÒÑÐγɣ¨42-72£¬13-42£©ºÍ£¨13-42£¬43-72£©

+ l+ q8 Z2 m& p8 d# k5 o%ÎÞ¹¦²¿·Ö

( ^' V: L- J( V% Hfor i=1:30

" C& o5 C3 N* D8 s4 j7 T
for j=1:30

* i+ K4 [* X8 L9 |' g0 A- M- u' ^
for k=j:30

5 I# {* n: P4 [# d- N( o+ x$ p( U" oif (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

" j( o8 B5 W' P0 x* z8 |& U: d5 I
elseif (j==k&&i==j)

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

                    temp=0;                   %thth

# T& L) r! `! W' x, s' E$ _for l=1:30

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

! A" e' |3 y4 t: M% d
end

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

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

+ O9 `8 W- @( U- I& 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# ~, R1 ^" Q% L- felseif (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);

4 P6 \5 A/ o9 v8 U% O: u
end

' R% z2 P# v$ p! [( e) F; O; y
end

) i9 ?6 F. q( ~6 b) rend

* x4 U4 W# M9 g, H) V# Fend+ d( O# h, R! z% ]2 T- E" A5 C) ?
%ÖÁ´ËÒÑÐγɣ¨13-42£¬13-42£©ºÍ£¨42-72£¬43-72£©

; w* X2 z' v5 jfor i=1:30

9 Y. N$ A& c. H( s% v. @# _8 Tfor j=1:30

4 y) R% X' s9 w8 ?1 e
for k=1:30

- K; p) d8 Y- p1 y2 i2 M& o0 u/ E
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

  D  X. C1 q4 i! [# p* A
elseif (j==k&&i==j)

                    temp=0;                   %thV

4 G/ ~6 i* n1 l+ efor l=1:30

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

9 v7 f6 @2 w. j* L8 G& j& q2 \0 o- Lend

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

2 |* `* \  F4 ~
elseif (j==i)

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

. H& {* N4 f5 E- N$ kelseif (k==i)

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

  Y, x, ~6 b- @# P: Z
end

$ s3 @) E5 x0 g% hend

+ V+ _& T$ J- S0 K
end

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

3 i  ~  o3 J" h$ P( |* hend9 K) @# o" x  n& Z4 t1 ^" H
%ÖÁ´ËÒÑÐγɣ¨42-72£¬13-42£©ºÍ£¨13-42£¬43-72£©

2 n: g1 d: |- J; \: M8 j8 Q+ v+ n
%HhÐγÉÍê±Ï

! w' i4 v& d8 k9 k%Fourth Step: Jg, Hg

    Jg=eye(42,42);

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

    Hg=zeros(72);

# ~* x" ?( n$ g' u- W) x: ~
%Calculation Jacobian&Hessian matrix END

0 d( ]. q  R" e7 K, s3 g; o
%%

4 E$ o/ I6 N- t* k. L) T' T
%Calculate Newton Iteration Îó²îµü´úÁ¿

) l' _. w6 r2 x9 F+ m& C- V8 c  n%Cal LX0-------------------------1

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

7 Q- M, u) T% c# T. b- f, m! ~%Cal LLam-------------------------2

    LLam0=h;

    pferr=max(LLam0);

  f6 |( N$ u0 `8 w; s5 f, G7 K3 D
%Cal LMU_MIN-------------------------3

    LMU_MIN0=g-sl-gmin;

% G2 C5 o( W* V7 |
%Cal LMU_MAX-------------------------4

    LMU_MAX0=g+su-gmax;

. ~  a/ K% d+ N+ f
%Cal Lsl-------------------------5

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

+ E* g6 d& Y7 V" y8 i0 h, Y2 \! w7 r
%CAl Lsu-------------------------6

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

5 W+ V" }% m# H; E0 V
%Calculate Newton Iteration Îó²îµü´úÁ¿ END!!!

* G9 P  v: c" S% O  n  ~
%%

1 |1 M# w9 _& D: R. V/ d, H
%Calculate Newton Iteration ·ÂÉäÐÞÕýÁ¿

+ V( |! L& i$ n3 D8 i: {( E%1st Step: ·ÂÉäÐÞÕýÁ¿delXaf,·ÂÉäÐÞÕýÁ¿delLamaf

) M8 C  a  @2 |& t: @# F: P%S1:H

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

6 _5 w% W: ]) i) c0 ?
for i=1:60

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

8 a* W5 D' D( K4 ~4 _* ?end

     tempgMUg1=0;

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

% M( @* \/ ^5 v0 o& w; N9 F1 C0 lfor i=1:42

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

5 {2 U+ _) z0 w0 D1 x% s
end

     tempgMUg2=0;

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

, _, n4 K( w" h+ i0 r7 F4 ~for i=1:42

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

" ?- R& Z% G# V8 h* q
end

     H=Hf-temp+tempgMUg1+tempgMUg2;

3 {) b6 a! ]/ x& V) M5 W! ?* |
%S2:·ÂÉäKESAaf

     tempgMUg1=0;

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

; H0 G+ M. F: o8 S* C, @  mfor i=1:42

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

+ ]$ w! n4 w2 V
end

     tempgMUg1=Jg*tempgMUg1;

     tempgMUg2=0;

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

& [3 c% @9 L3 @! t- G8 jfor i=1:42

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

1 z; F% k9 t4 p5 \/ s& Z, z+ {
end

     tempgMUg2=Jg*tempgMUg2;

     KESAaf=LX0+tempgMUg1-tempgMUg2;

1 N, [: U3 e0 ]" S) Q4 R%S3:JACOBIAN

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

  }4 d& r% `8 x+ Q0 _5 N" Y
%S4£ºRESULT

     RESULT=-[KESAaf;LLam0];

7 x, T0 f: o7 D0 N: w%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

% L& z# c7 E% O+ |( a; K

6 @; c4 Y% ?4 |- e%2nd Step: ·ÅÉäËɳڱäÁ¿ÐÞÕýÁ¿delslaf, delsuaf

* F; a  y9 Y6 X' q1 A
%S1: delslaf

     delslaf=Jg'*delXaf+LMU_MIN0;

4 d/ y% @0 `% ?- O
%S2: delsuaf

     delsuaf=-Jg'*delXaf-LMU_MAX0;

0 E( n# l8 f- t- Z* Z

9 J/ D1 o: a& ?# J. ~! J7 C, _
%3rd Step: ·ÅÉäÀ ¸ñÀÊÈÕ³Ë×ÓÐÞÕýÁ¿delMU_MINaf,delMU_MAXaf

" ]$ A# q6 i& B%S1: delMU_MINaf

     temp=0;

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

6 y- q6 k, ~8 a! [2 g% B8 T: F. \for i=1:42

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

9 j' W+ w' t3 K1 k3 I: ]end

     temp(13)=0;

     delMU_MINaf=temp;

& j" t+ a/ s. G5 @/ L* F%S2: delMU_MAXaf

     temp=0;

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

7 U7 f4 R: _$ D9 d3 K) xfor i=1:42

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

, S% Y- H% g4 S0 R& u1 _! j
end

     temp(13)=0;

     delMU_MAXaf=temp;

9 D$ ]% n* }; i

" T; b2 i+ u! Z* G; z& t. B
%Calculate Newton Iteration ·ÂÉäÐÞÕýÁ¿ END!!!!!!!!!!!!!!

1 s6 X& l2 l9 y! L/ A9 E

" v8 t* Y. g0 W: i%%

, Y( T* i) K+ J, i# ^* O%¼ÆËã·ÂÉäÔ Ê¼²½³¤ºÍ¶Ôż²½³¤ STEPpaf£¬ STEPdaf

; @0 X/ X4 ?, _9 I%S1: STEPpaf

. ]- j$ Q& p1 v# p$ L+ `
for i=1:42

4 t. S" X4 D# \
if (delslaf(i)~=0)

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

4 s& l: J  Y' b' q  welse

            temp1=Inf;

- r9 H% b( U, q+ f  y: g& v/ s
end

4 n' p! S1 Z# E7 m2 zif (delsuaf(i)~=0)

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

" }# f* I: k4 T1 R; J" Z: ^; Lelse

            temp2=Inf;

3 t% P2 C; c& W& C8 F1 g) N) W
end

4 V/ Q8 p, E; o2 k5 z0 x3 h! w
if temp1<temp2

            min=temp1;

5 w9 ~  U3 i" m- |. h
else

            min=temp2;

9 g1 X0 V3 E# H
end

( R2 m" G, r+ k, D2 pif min>1

            min=1;

4 ^- Q4 M8 c- iend

, D' Z: t% ?+ B, G( F
end

    STEPpaf=0.9995*min;

( x5 X! s& q9 V0 q" M- N8 R! E%S2: STEPdaf

6 E8 u% J. x% o$ t6 ?+ [for i=1:42

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

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

8 M; A, @- Y% P5 r  H$ Q3 qif temp1<temp2

            min=temp1;

0 D! K9 |# `+ Y3 M% U" ^else

            min=temp2;

; R2 m4 B, h# _' N* T8 a, ^- [6 oend

' ?8 Q( V- A9 S; n
if (i==13)

            min=0;

+ c9 C  G0 P% {7 S  \- uend

: ~5 W/ {9 V0 ^$ y2 d  w# Z
if min>1

            min=1;

2 t# j+ X0 W, [
end

: R/ M/ x( F1 z0 b, w, [0 Q
end

    STEPdaf=0.9995*min;

: J8 {; \( I% C+ T, T7 N
%&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!!!!!

5 T( w' _' r/ i7 l) j* v$ C" ?8 H. `

* ?) l' |1 V" {%&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 l/ D* I/ U5 K, B
if (ROUaf/ROU)^2<0.2

       SIGMA=(ROUaf/ROU)^2;

* G% ^: q) A( C8 E; Ielse

       SIGMA=0.2;

/ j0 x' v( r3 O* O
end

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

- v  Z# `- r9 S% {: \8 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;

8 p$ Z, e; s/ @) V%Calculate Newton Iteration &ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;

. o0 f9 w( B* B, ^%1st Step: &ETH;&pound;&Otilde;&yacute;&ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;delX,&ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;delLam

& n) V- O* i5 V& g%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;

; d" x  X% u9 D8 @# afor i=1:60

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

5 j9 x% A2 \" a% _
end

     tempgMUg1=0;

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

; T9 w: U: y: u) {+ ^for i=1:42

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

- V& G, V, n3 z! {+ ?( t% cend

     tempgMUg2=0;

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

7 d& E2 _0 }+ afor i=1:42

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

& T0 y3 [# T: D) C: N/ e: U" iend

     H=Hf-temp+tempgMUg1+tempgMUg2;

+ `: ]$ d7 e7 m4 b  c%S2:&ETH;&pound;&Otilde;&yacute;KESA

     tempgMUg1=0;

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

4 e7 V, D  j/ P( K! T+ tfor i=1:42

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

/ k$ S! i! y$ ^; T& D
end

     tempgMUg1=Jg*tempgMUg1;

     tempgMUg2=0;

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

: ?5 `* d3 a" Hfor i=1:42

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

: `+ _1 d% Y9 X' V: Oend

     tempgMUg2=Jg*tempgMUg2;

     KESA=LX0+tempgMUg1-tempgMUg2;

6 R  C6 T; G" [%S3:JACOBIAN

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

  x' |& `* |+ T2 V; A. e
%S4&pound;&ordm;RESULT

     RESULT=-[KESA;LLam0];

/ A0 Q+ i9 T1 v8 H( h% A$ E- M7 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

  _9 d. \; E( J; ^

4 O$ ~" y! t- D3 p4 I3 ?/ Z4 p6 C. t%2nd Step: &ETH;&pound;&Otilde;&yacute;&Euml;&Eacute;&sup3;&Uacute;±&auml;&Aacute;&iquest;&ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest;delsl, delsu

8 W/ A. Q7 s$ u: j: o$ G2 p
%S1: delsl

     delsl=Jg'*delX+LMU_MIN0;

  Y! Z- @4 o7 l! a& b& N. {$ R. z: F
%S2: delsu

     delsu=-Jg'*delX-LMU_MAX0;

% e# k( |$ _# z, u$ q9 d

, F- A2 ?8 `8 B: l1 |+ l% E" E
%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

6 U! p* b1 j$ J# b7 i( a' K%S1: delMU_MIN

     temp=0;

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

) _% S" C0 s0 e3 x6 N) @$ m( ^5 @for i=1:42

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

# b7 E( p8 {& H1 k/ E. b
end

     temp(13)=0;

     delMU_MIN=temp;

$ x1 b4 v0 Q9 r& p1 v3 Z
%S2: delMU_MAX

     temp=0;

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

8 ]4 f: O0 E: [$ Q* Z
for i=1:42

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

9 ]0 P% k. q+ I9 A) Jend

     temp(13)=0;

     delMU_MAX=temp;

2 h7 X" m% E7 F  C8 C5 R

7 a4 h: _$ g! U! Z6 `& f6 O( T% X%Calculate Newton Iteration &ETH;&THORN;&Otilde;&yacute;&Aacute;&iquest; END!!!!!!!!!!!!!!

/ N$ n8 K9 ]$ X+ z
%%

7 m1 X0 i) r, k# ?! i0 R4 A%&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

, U% K& V  @3 U" o( t! c%S1: STEPp

+ z% t. a8 r" B4 Z# z( Z8 Q3 ^: ~
for i=1:42

7 c9 l2 ^) N( |3 Lif (delslaf(i)~=0)

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

$ T# x9 Z# w9 aelse

            temp1=Inf;

8 b" v. h: B# k! r
end

! _* m- R2 c+ tif (delsuaf(i)~=0)

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

% }0 d6 A7 |5 l  c  Xelse

            temp2=Inf;

& O. C& I* G8 J
end

# |$ a' h: {8 c/ M
if temp1<temp2

            min=temp1;

% @& s' A  z3 C+ A" K. _: yelse

            min=temp2;

1 g0 Y1 X0 |# r9 z& v& x. J- }2 u, C, Fend

% r6 k) J, y! T; D  l3 R1 H
if min>1

            min=1;

4 w8 ~, h2 F7 m7 y
end

5 H0 D# i0 P% Y& l; g1 A/ A# {
end

    STEPp=0.9995*min;

! i* q* O7 C! }%S2: STEPd

0 |! F1 H7 [! V
for i=1:42

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

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

' j9 d1 g9 z! N+ O# m
if temp1<temp2

            min=temp1;

1 Q( Z. t# J) p% I
else

            min=temp2;

; \1 x. N6 d) C/ C) s: C" ~end

; C! l+ S( s' ^( _if (i==13)

            min=0;

8 \9 Y$ W4 P& R' Mend

7 `9 p% M8 _0 L$ f  T/ R# ]if min>1

            min=1;

  f% f0 L- Z( |/ i2 @8 Q: @- E7 w
end

1 G- C8 `$ v( [+ G( ~( Y0 x* Vend

    STEPd=0.9995*min;

& d0 w2 ?; x0 F+ ?8 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);

3 J; ?; Y) J9 n% Q%&cedil;ü&ETH;&Acirc;&Ocirc; &Ecirc;&frac14;±&auml;&Aacute;&iquest;&ordm;&Iacute;&para;&Ocirc;&Aring;&frac14;±&auml;&Aacute;&iquest; END!!!

6 w. {/ v" e7 E

/ |7 z+ I3 w# X, i
%&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));

. r" e. j! K' u1 Y* i1 W
%%

; ?! X4 |' a, u: F$ B( o%&Aring;&ETH;&para;&Iuml;&Ecirc;&Ccedil;·&ntilde;&sup3;&not;&sup1;&yacute;&micro;ü&acute;ú&acute;&Icirc;&Ecirc;&yacute;

    ik=ik+1;

  m3 k  \& k) A7 f! N8 e- ]# o
if ik>IterNumMAX

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

5 O- G# @/ D  h) J
break;

+ U, C/ f. L9 F7 Rend

/ K; k4 b4 J7 B. Y7 O

end

et = etime(clock, t1);

%%

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

* K6 n1 h+ ^9 m3 R% Q( C4 H

if (ik<=IterNumMAX)

9 K/ Y* Z8 l6 `5 E
%%

7 G6 y' F) S9 R# `# Z8 O" Y6 L+ H! l
%&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

' [. h+ ]2 h( ^/ m* o3 j, q: |& a* `

max=-10;min=10;

for i=1:ik

    temp=LLam0(i);

* \' j# j* g/ G& ~if (temp>max) max=temp;

7 c; A6 F# q% R. o: y+ V$ uend

$ I( a: |8 F$ z. ~8 c
if (temp<min) min=temp;

, d: Z" ]% m3 }0 H! z: Q0 w4 rend

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

  l: f6 {. v9 t. B0 kfor i = 1:30

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

0 l; r* `3 |9 `8 q' Fif (i<=6)

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

2 N! O3 z2 @; m7 F6 o+ Y4 i9 q" m
else

            fprintf('       -         -  ');

. z& n0 c" p5 m( x2 ?0 send

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

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

, L6 o% G( {) M  E" ~+ ]& ]
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');

- h. [! x8 N- x- w5 {" n
end

( A" b8 A  _+ |, K# d
end

* z7 w) d7 }4 ?6 b
%%

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