問題描述

現代電力系統分析課程留了用計算機進行潮流計算的作業，可以使用各種軟件實現，matlab就包含在其中，按老師的意思應該是使用matlab中的某個app。

可是我頭鐵，app是什么，當然要自己寫啊！

以此題為例

題目來自于 電力系統分析教材 華中科技大學版

整體思路

根據Y陣，PQ節點的P、Q值，PV節點的P、V值，平衡節點的δ、V值，利用極坐標的牛頓法，求出每個節點的δ、V值和P、Q值，和ΔP、ΔQ值，根據ΔP、ΔQ值判斷是否收斂，同時得出下一個P、Q值，反復迭代求解。

原理及步驟

輸入Y陣Ymat，PQ、PV節點參數Ps、Qs、Ps、V，平衡節點參數，此為目標值。

給出初始V、δ值，其中V節點、平衡節點為1中的值，在迭代過程中不再變化。

根據δ、U值求其對應的P、Q值。

根據P、Q、δ、V、Ymat求Jacobi矩陣J。

將復數矩陣Ymat進行分離，實部組成G陣，虛部組成B陣。根據如下公式計算J陣。

根據Ps與P，Qs與Q的差求出修正量DP、DQ。

根據J陣與DP、DQ，利用牛頓法公式

其中:

求Δδ、ΔV。

根據Δδ、ΔV、δ、V求得下一個δ、V。

重復3~8。

程序代碼

/chaoliujtz.m

%季天澤

%華中科技大學電氣與電子工程學院

%強電磁與新技術國家重點實驗室

clear

clc

% 節點必須按照PQ、PV、平衡節點來排序

% PQ節點有m個，共有n+1個節點，即PV節點有n-m個

% 每個PQ節點有1個關于Q的方程式，每個PV或PQ節點有1個關于P的方程式，平衡節點沒有方程式

% 所以P有n行，Q有m行，U與δ有n+1行

Ymat=[

1.042093-i*8.242876 -0.588235+i*2.352941 i*3.666667 -0.453858+i*1.891074

-0.588235+i*2.352941 1.069005-i*4.727377 0 -0.480769+i*2.403846

i*3.666667 0 -i*3.333333 0

-0.453858+i*1.891074 -0.480769+i*2.403846 0 0.934627-i*4.261590

];

z=max(size(Ymat));

%Ps、Qs代表實際的P、Q

Ps(1)=-0.3;

Ps(2)=-0.55;

Ps(3)=0.5;

Qs(1)=-0.18;

Qs(2)=-0.13;

n=size(Ps);

m=size(Qs);

n=n(2);

m=m(2);

%delta=[ones(1,10)*-1,0]

delta(z)=0;

%輸入U時，注意PV節點和平衡節點直接輸入的是目標值的標幺值，而PQ節點輸入的僅是迭代初值

U=[1,1,1.1,1.05];

output=[];

error=[];

G=real(Ymat);

B=imag(Ymat);

nn=10%將迭代次數設為10

for ii=1:nn

p1=(U*(G.*cos(delta-delta')+B.*sin(delta-delta')));

q1=(U*(G.*sin(delta-delta')-B.*cos(delta-delta')));

P=U(1:n).*p1(1:n);

Q=U(1:m).*q1(1:m);

clear p1 q1;

DP=Ps-P;

DQ=Qs-Q;

J = Jacobi(P,Q,Ymat,U,delta);

invJ=inv(J);

%修正方程

M=-invJ*[DP';DQ'];

Ddelta=M(1:n)';

DU=(diag(U(1:m))*M(n+1:m+n))';

Ddelta(n+1)=0;

DU(n+1)=0;

delta=delta+Ddelta;

U=U+DU;

error=[error;[DP,DQ,DU]];

output=[output;U,delta];

end

plot(1:nn,output)

/Jacobi.m

%季天澤

%華中科技大學電氣與電子工程學院

%強電磁與新技術國家重點實驗室

%

function J = Jacobi(P,Q,Y,U,delta)

% U,Y,S都是復矩陣，獲取之后進行拆分,m為PQ節點的個數，n為節點總數，n-1-m為PV節點數

% 節點按照PQ、PV、平衡節點來排序

% PQ節點有m個，共有n+1個節點，即PV節點有n-m個

% 每個PQ節點有1個關于Q的方程式，每個PV或PQ節點有1個關于P的方程式，平衡節點沒有方程式

% 所以P有n行，Q有m行

n=length(P);

m=length(Q);

G=real(Y);

B=imag(Y);

Q(n+1)=0;

P(n+1)=0;

for i=1:n

for j=1:n

if i~=j

H(i,j)=-U(i)*U(j)*(G(i,j)*sin(delta(i)-delta(j))-B(i,j)*cos(delta(i)-delta(j)));

elseif i==j

H(i,j)=U(i)*U(j)*B(i,j)+Q(i);

end

end

end

for i=1:n

for j=1:m

if i~=j

N(i,j)=-U(i)*U(j)*(G(i,j)*cos(delta(i)-delta(j))+B(i,j)*sin(delta(i)-delta(j)));

elseif i==j

N(i,j)=-U(i)*U(j)*G(i,j)-P(i);

end

end

end

for i=1:m

for j=1:n

if i~=j

K(i,j)=U(i)*U(j)*(G(i,j)*cos(delta(i)-delta(j))+B(i,j)*sin(delta(i)-delta(j)));

elseif i==j

K(i,j)=U(i)*U(j)*G(i,j)-P(i);

end

end

end

for i=1:m

for j=1:m

if i~=j

L(i,j)=-U(i)*U(j)*(G(i,j)*sin(delta(i)-delta(j))-B(i,j)*cos(delta(i)-delta(j)));

elseif i==j

L(i,j)=U(i)*U(j)*B(i,j)-Q(i);

end

end

end

J=[

H, N;

K, L

];

程序運行結果

數據迭代結果

數據迭代結果折線圖

結論

可以看出，用牛頓法求解電力系統潮流方程，僅需3步即可收斂，收斂速度極快。

你可以隨意更改我的代碼的任何部分，如果我的思路和代碼對你有所幫助，請點贊或評論，這對我非常重要，謝謝

總結

以上是生活随笔為你收集整理的牛顿法的matlab实现,利用matlab基础语句实现牛顿法matlab潮流计算的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。