matlab计算npv,哪位大神看看这个程序呀，参数下面给了

该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

K=0;

N=input('请输入节点数:N=');

NPQ =input('请输入节点数:NPQ=');

NPV =input('请输入节点数:NPV=');

Nl=input('请输入支路数:Nl=');

B1=input('请输入由各支路参数形成的矩阵:B1=');%[1 2 0.1+0.4i 0.01528i 1 1;1 3 0.3i 0 1.1 0;1 4 0.12+0.5i 0.0192i 1 0;2 4 0.08+0.4i 0.01413i 1 0]格式为某支路的首端号P，某支路末端号Q，且P

Y=zeros(N);

for i=1:Nl

if B1(i,6)==0

p=B1(i,1);q=B1(i,2);

else p=B1(i,2);q=B1(i,1);

end

Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));

Y(q,p)=Y(p,q);

Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;

Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;

end%求导纳矩阵

disp('导纳矩阵Y=');

disp(Y);G=real(Y);B=imag(Y);

Kmax=input('\n\n 请输入最大迭代次数后回车(可从零开始) Kmax=\n');%防止发散的情况可以及时跳出

small=10^(-4); %ε不能太小

Pnode=[ 0.2 -0.45 -0.4 -0.6 ];

Qnode=[ 0.2 -0.15 -0.05 -0.1 ];

Vnode=[ 0 0 0 0 ];% 书上给定的有功无功电压初值

e=[ 1.06 1.0 1.0 1.0 1.0 ];

f=[ 0 0 0 0 0 ];% 正常给定的初值

for K=0:Kmax,

%计算△W。△W=【△P1,△Q1,△P2,△Q2,……,△P,△V^2(平方的意思)】(平衡节点不参与迭代)

%用书上134页及之后计算NPQ个PQ节点的△P、△Q。算法：先初始化，再不断地累减

%初始化,得△P=【Pnode(1),Pnode(2),……，Pnode(NPQ)】PQPV节点个数上面要输入

for i=1:NPQ,

dP(i)=Pnode(i);

dQ(i)=Qnode(i);

end

%不断地累减,得△P和△Q ，就是公式啦

for i=1:NPQ,

for j=1:N,

dP(i)=dP(i)-e(i)*G(i,j)*e(j)+e(i)*B(i,j)*f(j)-f(i)*G(i,j)*f(i)-f(i)*B(i,j)*e(j);%△P(4-86)

dQ(i)=dQ(i)-f(i)*G(i,j)*e(j)+f(i)*B(i,j)*f(j)+e(i)*G(i,j)*f(i)+e(i)*B(i,j)*e(j);%△Q一样

end

end

%计算NPV个PV节点的△P、(△V)^2 (其中△P的计算同上)算法：先初始化，再不断地累减，电压就继续公式

%初始化

for i=(NPQ+1):(N-1),

dP(i)=Pnode(i);

end

%不断地累减,得△P和(△V)^2

for i=(NPQ+1):(N-1),

for j=1:N,

dP(i)=dP(i)-e(i)*G(i,j)*e(j)+e(i)*B(i,j)*f(j)-f(i)*G(i,j)*f(i)-f(i)*B(i,j)*e(j);

dV2(i)=Vnode(i)^2-e(i)^2-f(i)^2;

end

end

%都是公式，不用太在意，PQ，PV节点分开算的(P134)

%组合成△W=【△P1,△Q1,△P2,△Q2,……,△P,△V^2】(平衡节点不参与迭代)

%将△P间隔地赋入△W

a=1;

for i=1:(N-1),

dW(a)=dP(i);

a=a+2;

end

%将△Q间隔地赋入△W

a=2;

for i=1:NPQ,

dW(a)=dQ(i);

a=a+2;

end

%将△V^2间隔地赋入△W

a=NPQ*2+2;

for i=(NPQ+1):(NPQ+NPV),

dW(a)=dV2(i);

a=a+2;

end

%判断是否小于ε，注意判断方法是△P，△Q的绝对值足够小

if max(dW)

fprintf('\n 迭代是收敛的。第%d次迭代后终止迭代。\n',K-1);

break;

end

%计算雅克比矩阵各元素(书上135页)。算法：先初始化，再不断地累减

J=zeros(N-1);

%初始化

for i=1:N-1,

for j=1:N-1,

if i<=NPQ %求雅克比矩阵中PQ节点的元素

if i==j%i=j表示对角元素

J(2*i-1,2*j-1) = -G(i,i)*e(i)-B(i,i)*f(i);

J(2*i-1,2*j ) = B(i,i)*e(i)-G(i,i)*f(i);

J(2*i ,2*j-1) = B(i,i)*e(i)-G(i,i)*f(i);

J(2*i ,2*j ) = G(i,i)*e(i)+B(i,i)*f(i);

elseif i~=j%非对角

J(2*i-1,2*j-1) = -G(i,j)*e(i)-B(i,j)*f(i);

J(2*i-1,2*j ) = B(i,j)*e(i)-G(i,j)*f(i);

J(2*i ,2*j-1) = B(i,j)*e(i)-G(i,j)*f(i);

J(2*i ,2*j ) = G(i,j)*e(i)+B(i,j)*f(i);

end

end

if i>NPQ %求雅克比矩阵中PV节点的元素

if i==j

J(2*i-1,2*j-1) = -G(i,i)*e(i)-B(i,i)*f(i);

J(2*i-1,2*j ) = B(i,i)*e(i)-G(i,i)*f(i);

J(2*i ,2*j-1) = -2*e(i);

J(2*i ,2*j ) = -2*f(i);

elseif i~=j

J(2*i-1,2*j-1) = -G(i,j)*e(i)-B(i,j)*f(i);

J(2*i-1,2*j ) = B(i,j)*e(i)-G(i,j)*f(i);

J(2*i ,2*j-1) = 0;

J(2*i ,2*j ) = 0;

end

end

end

end

%不断地累减

for i=1:NPQ,

for k=1:N,

J(2*i-1,2*i-1) = J(2*i-1,2*i-1)-G(i,k)*e(k)+B(i,k)*f(k);

J(2*i-1,2*i ) = J(2*i-1,2*i )-G(i,k)*f(k)-B(i,k)*e(k);

J(2*i ,2*i-1) = J(2*i ,2*i-1)+G(i,k)*f(k)+B(i,k)*e(k);

J(2*i ,2*i ) = J(2*i ,2*i )-G(i,k)*e(k)+B(i,k)*f(k);

end

end

%根据△W=-J*△V, 得△V书135页

dV=(-J)\dW';

%用dV对e、f进行修正,并得到复数表示的V

for i=1:(N-1),

e(i)=e(i)+dV(2*i-1);

f(i)=f(i)+dV(2*i );

V(i)=complex(e(i),f(i));

end

V(N)=complex(e(N),f(N));

format long g;

fprintf('\n 第%d次迭代(共%d个独立节点)',K,N);

V

end

N=5

NPQ=4

NPV=1

NL=7

B1=[1 2 0.02+0.06i 0 1 0;1 3 0.08+0.24i 0 1 0;2 3 0.06+0.18i 0 1 1;2 4 0.06+0.18i 0 1 1;3 4 0.01+0.03i 0 1 0;4 5 0.08+0.24i 0 1 0;2 5 0.04+0.12i 0 1 1;]

Kmax=50