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

莫言wan · 发表于 2013-7-7 15:10:59

马上加入，结交更多好友，共享更多资料，让你轻松玩转电力研学社区！

您需要登录才可以下载或查看，没有账号？立即加入

×

如何修改matpower的目标函数或者重建一个目标函数重新运行一下？求大神指教！！！

redplum · 发表于 2013-7-7 15:52:14

有没有例程贴上来看一下

qcsun · 发表于 2013-7-9 10:11:13

你要把它的程序读懂，相应部分都要修改的

莫言wan · 发表于 2013-7-10 10:32:27

回复 2# redplum

首先谢谢指教
下面是这个程序，指教一下那里需要修改啊

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)

: \6 c1 y% J- d- z
%%

6 ^1 h% e/ z6 D/ b2 x$ |7 k E
%Calcute h,g matrix

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

errArr=[errArr;ROU;];

SIGMA=0;

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

, n v3 K* X8 m( V, `
for i=1:30

temp=0;

! e) _- i$ z7 w8 Xfor j=1:30

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

: J1 q1 m% X! m7 k. `0 p
end

. S, g, B/ U" q7 L% a" v3 P
if (i>6)

tPg=0;

. a1 F& O9 i+ Q3 F' j3 [- H4 }else

tPg=Pg(i);

0 M+ y& `. V. W3 aend

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

0 x: z6 K2 B) w$ B3 e' ^end

8 L1 z6 }1 M2 e$ b
for i=1:30

temp=0;

* B4 g; N" H6 K0 m8 d7 X8 i3 k8 ufor j=1:30

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

) z5 h$ J# ~9 }2 qend

; \8 k! a5 J) U$ Z ]2 u
if (i>6)

tQg=0;

3 v, {; L5 |% Y# Y7 a
else

tQg=Qg(i);

. L% e* w# P/ {% V( Y. g
end

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

; g$ i5 S7 E( r3 e/ c! P0 wend
N# ~) V$ G6 d s" l$ Y% Cal h END

% p9 n% q+ d+ o! J9 ^5 i, x

% M4 G8 Q$ r+ U1 u6 ? b
for i=1:6

g(i)=Pg(i);

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

& Z4 o8 U$ D, ?5 O+ Z
end

; I: ~+ N: U- n3 q: D$ S
for i=1:30

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

& Z% J9 ^' n' i! cend! @* ~3 t4 f [, N+ D" U; ]( c3 ^
% Cal g END

1 ?7 G& j( ]6 `) z0 d$ w# m
%Calcute h,g matrix END

4 E$ O0 K. z0 P' e%%

& i1 p/ H. k( Z) P%Calculate Jacobian&Hessian matix

5 @# q) w4 ?! D+ N3 o, j. l%First Step: Jf,Hf

" i% f3 r i6 E( N, X) O' ?( { pfor i=1:6

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

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

+ g, K1 c( g" o4 y) ]end

; f3 u! \+ E% i6 c" P$ W) I7 Y: |1 l%Second Step: Jh, hÎªµÈÊ½Ô¼Êø

) O" a: z! e, m/ o
for i=1:6 %Ç°6ÐÐ¶ÔPgÇóµ¼£¬ÓÉ´ËÒÑÇó³ö

Jh(i,i)=1;

2 N& t2 K6 D! [# B1 X4 U
end

8 h' @- x t* v+ e3 q+ \
for i=7:12 %7-12ÐÐ¶ÔQgÇóµ¼£¬ÓÉ´ËÒÑÇó³ö

Jh(i,i+24)=1;

9 F- \3 ]( _+ H) M
end

1 Y5 ]7 b4 z/ {. _4 s
for i=1:30 %ÐÎ³É13-42ÐÐµÄ1-60ÁÐ

) ~* N0 ^+ r6 Q
for j=1:30

tempVp=0;

tempVq=0;

: o) o/ V# s) M7 z/ ]! M& l* ~
if (j==i)

/ D, k1 J; Q, L+ r( L
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));

0 A5 g$ f2 j3 Z3 O4 _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));

; X: f: v1 a+ g, Gelse

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

% Z" @$ M- J* W: g. d
end

( c0 m8 p9 U! }. Vend

: R6 v8 Z2 t. L: z
end

1 j* h9 g# I& Z: j' _7 c
for i=1:30 %ÐÎ³É43-72ÐÐµÄ1-60ÁÐ

1 R% L# [1 O$ {, |9 F9 \( {
for j=1:30

tempVp=0;

tempVq=0;

9 `$ o0 Q, ]2 lif (j==i)

- Q( V( \, N+ t E4 j! n
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));

# u/ Y0 w; v2 i4 B
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;

9 w# g7 _7 e! L* w# |5 F) s
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));

- I2 w7 f6 B5 L4 T# x: z2 v" mend

( q9 O* z* N6 m/ V% z; q: q$ Qend

& J1 F t$ C1 @$ d
end

) o0 V A/ b5 I' {: E8 p8 A! W- x8 T%Third Step: Hh

6 G" V! O& o1 R) |
%ÓÐ¹¦²¿·Ö

* \, o/ |/ Q9 f/ J2 d/ ufor i=1:30

8 ~1 l, u6 n8 l
for j=1:30

% l: H4 `- K1 I" @
for k=j:30

) C4 `$ b0 ~1 }# Q2 Q [4 n
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

9 `0 A' C% Z7 Q2 V& ~elseif (j==k&&i==j)

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

temp=0; %thth

$ y, P8 c* Q- l& s qfor l=1:30

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

) X( N$ v5 D( d8 }: r4 j1 L/ R
end

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

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

. u! Q$ J- z5 m N2 j, u
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);

' o$ V- k0 E+ @( A0 N2 ^
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);

( ^1 ]9 Z, T) l/ }2 n4 j
end

( |1 H- _* Z4 t, uend

; V- O7 W6 Z: B: v$ K9 p: b- B1 \end

* P7 q% ]5 R5 F3 B$ R: v6 D; Dend# {6 n0 D5 ^+ k: {+ W
%ÖÁ´ËÒÑÐÎ³É£¨13-42£¬13-42£©ºÍ£¨42-72£¬43-72£©

2 b) h5 I& ? jfor i=1:30

( _# E& Z' j5 m
for j=1:30

5 J) R7 M& C7 e: u5 ^4 W4 rfor k=1:30

4 b: e. G3 J0 A" H9 D; G( [# ~
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

/ S- Q* H$ a: Z: C1 L4 W
elseif (j==k&&i==j)

temp=0; %thV

$ g1 i4 C* v! ?2 s- h* z; s
for l=1:30

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

* k: N3 i& h1 _$ B# q
end

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

- }- J+ e* L; d7 c" c
elseif (j==i)

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

# x: x% X8 c) n7 D2 f5 B
elseif (k==i)

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

! a7 T2 t; B' a5 Q# |7 tend

- W$ h5 j) B4 Mend

% g! n K8 ]1 T2 k% }" oend

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

4 e3 u" e, V/ I# Bend
- H/ I! x4 Z- ]6 S%ÖÁ´ËÒÑÐÎ³É£¨42-72£¬13-42£©ºÍ£¨13-42£¬43-72£©

1 Z: c0 g3 j3 q%ÎÞ¹¦²¿·Ö

5 Q" Z7 N. n" E$ Z8 {
for i=1:30

) E' S0 a+ w4 m& }2 Rfor j=1:30

; S2 X% r: M3 t- I
for k=j:30

( r9 x, O, @" P" Q
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

& ~ x1 t& X; T; g$ D( x9 `/ l
elseif (j==k&&i==j)

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

temp=0; %thth

1 I" ~. L; w) o- A9 h
for l=1:30

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

3 g/ }9 J2 H8 A( I8 W: Mend

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

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

+ R. I0 X/ x, Relseif (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);

) ]6 I7 n c: t! P# Q. }- v
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);

$ Q1 A" j; J" @; H1 e
end

$ @* ^% g, Z7 f6 j- a6 {
end

8 o2 x# G" Z$ |( W% I* S! M
end

4 _. d" Y' f, @" k
end
. z" D+ J7 c( Y2 l! _ h9 k%ÖÁ´ËÒÑÐÎ³É£¨13-42£¬13-42£©ºÍ£¨42-72£¬43-72£©

# _) V- |& s) i/ a* u* o
for i=1:30

9 r3 c: M" l( B5 W5 I/ `, qfor j=1:30

! p& q$ Q% J+ `5 |3 u
for k=1:30

* @6 H/ ]% s6 A% j# sif (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

4 `3 ~9 Q! R9 {' J+ p! P! ]
elseif (j==k&&i==j)

temp=0; %thV

1 I$ z- Z8 G0 M) \% F* ^
for l=1:30

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

8 y$ H; N! C; _% N8 b5 oend

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

: ^" z0 y6 C" r. m% N
elseif (j==i)

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

* r8 n' C% p @" o: Delseif (k==i)

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

0 p- v* G8 P% N0 t
end

/ ~+ z1 J) W( q3 _9 f7 u
end

7 K+ r" }9 j5 O, O) ~end

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

9 O& F- E9 d( G$ ]3 T' |7 zend) Z/ @& _+ f% O+ j" T9 o* j- @' j
%ÖÁ´ËÒÑÐÎ³É£¨42-72£¬13-42£©ºÍ£¨13-42£¬43-72£©

3 ~2 ?$ `6 B9 E; C+ U' Q3 i" G6 i%HhÐÎ³ÉÍê±Ï

8 x5 c9 q% e) }/ J! {7 b%Fourth Step: Jg, Hg

Jg=eye(42,42);

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

Hg=zeros(72);

" y: J" z: N" w! e
%Calculation Jacobian&Hessian matrix END

. G" i$ J) O9 N% m( G
%%

/ i4 A0 Z8 t2 B8 N# W3 C%Calculate Newton Iteration Îó²îµü´úÁ¿

]1 v" h9 h: v' ~5 f& _2 [) [( }%Cal LX0-------------------------1

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

; y+ |- |7 t2 ?" b4 G: Q: |$ h$ q
%Cal LLam-------------------------2

LLam0=h;

pferr=max(LLam0);

0 v; ^" {0 @( |0 N
%Cal LMU_MIN-------------------------3

LMU_MIN0=g-sl-gmin;

0 l7 O% u% \% ^. w8 M8 q) _1 a%Cal LMU_MAX-------------------------4

LMU_MAX0=g+su-gmax;

. J9 ^" N( ~2 a4 Q* T%Cal Lsl-------------------------5

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

: p$ k! c. |( y) X% r
%CAl Lsu-------------------------6

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

) {0 H' f/ v, `- T0 o& v& T%Calculate Newton Iteration Îó²îµü´úÁ¿ END!!!

) m$ C! ?# z k- }; H" r%%

$ y# P2 P% f$ F( g1 K" a# N
%Calculate Newton Iteration ·ÂÉäÐÞÕýÁ¿

) z, M- R! O% j- z1 f- I- L%1st Step: ·ÂÉäÐÞÕýÁ¿delXaf,·ÂÉäÐÞÕýÁ¿delLamaf

% O+ I" p) z8 F; w
%S1:H

temp=0; %Ö®Ç°ÒÑÓÃ¹ýÁÙÊ±±äÁ¿temp£¬ÔÚ´ËÇåÁã

+ C; t3 \0 h( q& j1 e7 |
for i=1:60

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

/ B- n/ u* n# v* I4 F1 xend

tempgMUg1=0;

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

- P1 y' g. V% R+ c" Z' S
for i=1:42

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

J7 ]: m6 {+ B1 Fend

tempgMUg2=0;

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

/ p, y9 Q8 N/ @! C5 ]9 d( Hfor i=1:42

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

( ~9 R' n' p% Q. [' rend

H=Hf-temp+tempgMUg1+tempgMUg2;

+ i6 z0 }6 f) |! i2 e' W
%S2:·ÂÉäKESAaf

tempgMUg1=0;

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

- n+ I7 i f) B) b' ?; t
for i=1:42

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

5 c. C( F4 A! x; s- {1 Tend

tempgMUg1=Jg*tempgMUg1;

tempgMUg2=0;

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

# @& d2 r5 C* ufor i=1:42

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

2 F$ q. S) K0 y+ [, j2 ]
end

tempgMUg2=Jg*tempgMUg2;

KESAaf=LX0+tempgMUg1-tempgMUg2;

9 a! d% I: H( ^3 J%S3:JACOBIAN

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

0 z. q7 j3 |7 `2 o
%S4£ºRESULT

RESULT=-[KESAaf;LLam0];

: C0 R: b, R! p( C6 n D% J
%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

2 J4 ?, ^ }/ v% s4 d$ `+ E" n. {' v

: g( Q% y% ]0 k* S- \
%2nd Step: ·ÅÉäËÉ³Ú±äÁ¿ÐÞÕýÁ¿delslaf, delsuaf

# ]$ B' }, ~+ ]7 I* S2 {8 J9 ?: [; u
%S1: delslaf

delslaf=Jg'*delXaf+LMU_MIN0;

) e' X, P# y: Q* r: Y%S2: delsuaf

delsuaf=-Jg'*delXaf-LMU_MAX0;

i5 N9 y$ n" D9 W, ? j' F

0 ]5 U! s5 p1 F: U; u" U' ]) ]
%3rd Step: ·ÅÉäÀ ¸ñÀÊÈÕ³Ë×ÓÐÞÕýÁ¿delMU_MINaf,delMU_MAXaf

: c4 r' U# D/ x%S1: delMU_MINaf

temp=0;

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

4 {6 R, X3 W* E! E
for i=1:42

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

c) m' g" F/ i8 x3 }, t- p
end

temp(13)=0;

delMU_MINaf=temp;

0 y$ L+ i3 W3 D) B( S( Y) I& T%S2: delMU_MAXaf

temp=0;

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

, Y% t, R2 v) }# A4 P) G7 Bfor i=1:42

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

* G5 ~6 V6 b2 Z4 R2 ~" |* c
end

temp(13)=0;

delMU_MAXaf=temp;

1 }: e0 S, J; \" ?

7 Z# r( \! n( B%Calculate Newton Iteration ·ÂÉäÐÞÕýÁ¿ END!!!!!!!!!!!!!!

: b0 c7 y: R( J ~, \8 s

8 |7 g& [. o% v* e% |
%%

: c. m: R) F2 w$ [% }%¼ÆËã·ÂÉäÔ Ê¼²½³¤ºÍ¶ÔÅ¼²½³¤ STEPpaf£¬ STEPdaf

( ?; `- K4 V$ W0 B) }9 v4 D9 m%S1: STEPpaf

. b2 \2 s- O: U7 Hfor i=1:42

! f2 R+ w0 O6 c4 F: V
if (delslaf(i)~=0)

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

* a; `% h0 b. T8 m+ {, |else

temp1=Inf;

& j7 ^* T* P3 `, r( F& C- P9 c( E. rend

z8 h0 J; U3 l2 G1 ?" nif (delsuaf(i)~=0)

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

; p* Z( c& U! o o, V" y
else

temp2=Inf;

" d7 W4 u- E9 o: E5 dend

& n9 m. L1 X* l9 U& _$ T
if temp1<temp2

min=temp1;

. z7 R; {- ], C {else

min=temp2;

& Q* o7 Z: h4 Gend

& Q. Y6 Y6 x8 Y3 ^6 ^6 B; e0 c# oif min>1

min=1;

( s: {1 _) [9 \4 M& R5 Send

6 K* m/ b8 f" q; Qend

STEPpaf=0.9995*min;

% G% E1 j( {" A
%S2: STEPdaf

, l+ E* N' ` A7 n( `! S Yfor i=1:42

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

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

4 ?: W* ^6 G# E* O" R3 A+ oif temp1<temp2

min=temp1;

- f2 K! u. S9 X/ a4 ?6 m
else

min=temp2;

- e6 B0 w1 E: D0 R0 D+ ~- j% h
end

6 ]2 s" b+ E4 ^
if (i==13)

min=0;

j: h5 L4 Z! y! l. Bend

, e. N% N& f* T1 e9 E& e$ Lif min>1

min=1;

% d. ~9 E5 M* F! G& h, C$ q o+ [% M
end

, q) x/ U: s6 X1 Q/ o
end

STEPdaf=0.9995*min;

2 q! \; s* _) Q
%¼ÆËã·ÂÉäÔ Ê¼²½³¤ºÍ¶ÔÅ¼²½³¤ STEPp£¬ STEPd END!!!!!

$ D2 j' x& f3 w/ u' T6 X0 _

1 i0 d8 i- X: p" ^' F6 L%¼ÆËã·ÂÉä¶ÔÅ¼ÒòÊý¼°³Í·£Òò×Ó%

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

5 s' ?' r$ H; zif (ROUaf/ROU)^2<0.2

SIGMA=(ROUaf/ROU)^2;

# s% V/ L% K' s4 X2 Celse

SIGMA=0.2;

4 H+ v% v% q4 G! d* `end

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

2 w& \4 q" ?9 ^% \! t5 {- q0 X%¼ÆËãÍê±Ï%

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;

, I2 u Y5 ^3 O$ V* u9 V& A1 F
%Calculate Newton Iteration ÐÞÕýÁ¿

S& S: y2 u$ |5 n%1st Step: Ð£ÕýÐÞÕýÁ¿delX,ÐÞÕýÁ¿delLam

: {( L2 Y' U7 w; G8 x4 U4 \: ], H
%S1:H

temp=0; %Ö®Ç°ÒÑÓÃ¹ýÁÙÊ±±äÁ¿temp£¬ÔÚ´ËÇåÁã

& l0 m, I. X; Z1 f5 [- wfor i=1:60

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

+ W) F. N: R6 aend

tempgMUg1=0;

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

1 E, h& j W7 @& f. c, ~: s. B0 Dfor i=1:42

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

1 z7 O1 n# `+ dend

tempgMUg2=0;

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

% S5 j {, Z2 X& F
for i=1:42

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

- t3 Q/ b" \/ g( ~6 \8 yend

H=Hf-temp+tempgMUg1+tempgMUg2;

7 S+ Q! ~* G$ x8 Q4 o%S2:Ð£ÕýKESA

tempgMUg1=0;

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

( E" M w# ^& P- \8 A6 D
for i=1:42

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

& f, ^: l' @1 l8 W8 Q3 b
end

tempgMUg1=Jg*tempgMUg1;

tempgMUg2=0;

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

1 ~$ v' r- _0 n+ ?5 z$ qfor i=1:42

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

! S8 e, j5 |6 v" l8 o
end

tempgMUg2=Jg*tempgMUg2;

KESA=LX0+tempgMUg1-tempgMUg2;

' s# Q' `# C& m; O1 x
%S3:JACOBIAN

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

- [0 l. F6 v! h$ Z' Y
%S4£ºRESULT

RESULT=-[KESA;LLam0];

! _' Y+ Z! ~" [+ C8 Y" G%S5:Cal Ð£ÕýÐÞÕýÁ¿delX,ÐÞÕýÁ¿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

" W6 P5 |/ n0 E/ d+ @: f) g8 c3 P

' k T8 w, @& k0 B4 z* W" {( G$ R9 o
%2nd Step: Ð£ÕýËÉ³Ú±äÁ¿ÐÞÕýÁ¿delsl, delsu

& L8 U+ o) M5 e' ]%S1: delsl

delsl=Jg'*delX+LMU_MIN0;

4 i- g3 Y9 Q; m3 v
%S2: delsu

delsu=-Jg'*delX-LMU_MAX0;

" C& U0 Y9 B1 g, [

: W/ S5 d1 S B( `; G7 ?%3rd Step: Ð£ÕýÀ ¸ñÀÊÈÕ³Ë×ÓÐÞÕýÁ¿delMU_MIN,delMU_MAX

9 F" L$ h" h: ^/ W' N6 `1 ~' N
%S1: delMU_MIN

temp=0;

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

! N+ r/ m, l$ z+ _for i=1:42

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

1 ^# |( I$ g6 k9 e+ v: x) m
end

temp(13)=0;

delMU_MIN=temp;

" V6 p; G6 t$ k# C$ ~1 E+ Z; k%S2: delMU_MAX

temp=0;

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

3 p& t i8 x) G3 nfor i=1:42

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

/ B3 O4 ^% ?/ Nend

temp(13)=0;

delMU_MAX=temp;

2 [5 U( a5 b) @# X0 T

5 P% ^" \; X3 a%Calculate Newton Iteration ÐÞÕýÁ¿ END!!!!!!!!!!!!!!

6 C5 V% _# I0 I3 D%%

2 z: }) e' m* L. J8 M' n%¼ÆËãÐ£ÕýÔ Ê¼²½³¤ºÍ¶ÔÅ¼²½³¤ STEPp£¬ STEPd

, S* U7 ~4 Z' ]# c
%S1: STEPp

( u" m9 r& i2 D- w
for i=1:42

, O* D5 S! S& g5 V }+ h% m& T+ R, b
if (delslaf(i)~=0)

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

5 O4 o3 ~ u7 J4 c9 k: R
else

temp1=Inf;

! C/ j" x, k8 {& f$ Y5 K" b$ b, k
end

: A( }' }9 p) Z' x1 q; d0 s
if (delsuaf(i)~=0)

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

' G; T. n% U- V5 D) w5 }2 Z4 q# U
else

temp2=Inf;

; W5 f: v2 Z5 Q3 Mend

/ I' E4 R' D& Z" n+ }( j3 v* Tif temp1<temp2

min=temp1;

: d- U; F5 o5 d
else

min=temp2;

2 b/ L% a" O0 x( U2 \, ~5 R! |7 oend

( ^+ P. E0 d8 nif min>1

min=1;

/ A1 Q* R. F4 _: iend

4 ~- K$ ] B+ ^/ [4 C
end

STEPp=0.9995*min;

+ ` T% e4 @ k/ Y
%S2: STEPd

3 j9 \6 S2 k. i' I. S: afor i=1:42

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

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

) ~3 Y# d2 l. q6 s
if temp1<temp2

min=temp1;

3 b* d* W# k, z" f: x6 G6 \* ^
else

min=temp2;

. v" ?0 k$ r) j0 [6 h7 [) `( D
end

% ? A R# f3 L1 Bif (i==13)

min=0;

- {6 u2 x# G1 R/ k
end

5 V! W [, h7 _* E& V! S: `if min>1

min=1;

2 P2 P. w& T1 h4 i! s9 H* f8 i& ~end

4 W- ~7 x' S! N2 d8 O# j8 l" r
end

STEPd=0.9995*min;

; h* L( G. K0 S# q8 [5 r
%¸üÐÂÔ Ê¼±äÁ¿ºÍ¶ÔÅ¼±äÁ¿

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

, j* K) \, S$ T; M5 w/ @& x& ?$ ?%¸üÐÂÔ Ê¼±äÁ¿ºÍ¶ÔÅ¼±äÁ¿ END!!!

V, }* I% ~8 d% y

" w {9 e `( l7 t%¼ÆËã¶ÔÅ¼¼äÏ¶%

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

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

( M7 m" T8 B" s6 n
%%

4 r2 }7 ^3 [0 c0 w/ p# R. `%ÅÐ¶ÏÊÇ·ñ³¬¹ýµü´ú´ÎÊý

ik=ik+1;

+ L; i m1 F+ H/ k* L# s% `if ik>IterNumMAX

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

7 s6 N7 N- G# E } |& x) }- I+ b
break;

: J0 H. q4 M* T1 m( K
end

4 { N# q- J* V+ J

end

et = etime(clock, t1);

%%

%Êä³ö²¿·Ö

( {2 ?! I8 a. r1 b8 `

if (ik<=IterNumMAX)

) i% C, @3 ? d. |' h% T9 [. g%%

) F& m7 P8 Z! w4 C
%ÇóµÃ×î´ó×îÐ¡µÈÊ½Îó²î

F=0;

for i=1:6

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

end

3 R/ X! u7 R+ \! O# T

max=-10;min=10;

for i=1:ik

temp=LLam0(i);

" i, q' ]0 M) t4 I% r. E7 n1 @6 uif (temp>max) max=temp;

7 ?7 D" f8 r% C! b. a! j, I
end

3 q6 S1 C2 P3 H: R
if (temp<min) min=temp;

( e2 |. @4 O: b( ` ~end

end

max

ik

%%

%Êä³ö½á¹û

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

, A/ b6 v0 A ^/ H$ `2 mfor i = 1:30

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

2 E) u6 e: \7 M% X9 b& j& n
if (i<=6)

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

3 ?6 A% P# i2 p/ O" `
else

fprintf(' - - ');

- a% X0 Z# |/ _ L& h Q
end

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

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

( ~% ~4 e/ t- d b& Wif (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');

& k. g% E6 h- i* m, y" ^0 c
end

' ^2 U6 H) K- }$ j9 a, F6 Q. G+ t5 z
end

, Q$ |% q( ]8 q1 d5 s8 _
%%

%%Êä³öÇúÏßÄâºÏ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

莫言wan · 发表于 2013-7-10 10:33:53

回复 3# qcsun

谢谢，忙看一下哪里需要修改嘛，程序在二楼

mingyu · 发表于 2013-7-10 20:00:11

没看懂这个是做什么的。

当年依旧 · 发表于 2013-10-28 15:35:53

matpower里有自带的m文件，直接调用即可

账号		自动登录	找回密码
密码			立即加入

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

马上加入，结交更多好友，共享更多资料，让你轻松玩转电力研学社区！

相关帖子

浏览过的版块