matlab.txt

这是一段matlab程序

%---------------------------------------------------%
%   很早以前编写的小脚本，用MATLAB运行即可。        %
%   希望大家指教，可以到我的主页上留言。            %
%   http://hi.baidu.com/lyz_4534                    %
%---------------------------------------------------%





function lyz1()
%===================原始数据录入模块=====================%
error=1e-5;
max_n=10000;
balance_power_node=1.06+0i;
power=[ 0 , 0.2+0.2i , -0.45-0.15i, -0.4-0.05i , -0.6-0.1i];

%以下是两组数据都可以达到收敛,其中一个是书上(东南大学-稳态分析164页 4-3例题， 

%这两组数据分别为导纳或阻抗的形式，本程序可以自动识别。

put_in1=[%东南大学-稳态分析164页的4-3例题参数（导纳形式）
               1,3,                 1.250- 3.750i;
               3,4,                 10.00-30.000i;
               4,5,                  1.25- 3.750i;
               2,5,                  2.50- 7.500i;
               2,4,                 1.667- 5.000i;
               2,3,                 1.667- 5.000i;
               1,2,                 5.000-15.000i];
put_in2=[%一组我自己附加的参数（阻抗形式），供参考验证
               1,3,                  0.040+0.120i;
               3,4,                  0.005+0.015i;
               4,5,                  0.040+0.120i;
               2,5,                  0.020+0.060i;
               2,4,                  0.030+0.090i;
               2,3,                  0.030+0.090i;
               1,2,                  0.010+0.030i];
put_in=put_in1;       %选择其中的4-3例题参数进行计算
%=================节点导纳矩阵形成模块===================%
mn=size(put_in);
sum_node=max(max(put_in(:,1)),max(put_in(:,2)));
if             sum(imag(put_in(:,3)))>=0;
               put_in(:,3)=(linspace(1,1,mn(1))')./put_in(:,3);
end;
admittance=zeros(sum_node,sum_node);
for i=1:mn(1);
               if put_in(i,1)> put_in(i,2);
                   temp=put_in(i,1);
                   put_in(i,1)=put_in(i,2);
                   put_in(i,2)=temp;
               end;
end;
for node=1:sum_node;
               y_temp=0;
               for i=1:mn(1);
                   if put_in(i,1)==node||put_in(i,2)==node;
                       y_temp=y_temp+put_in(i,3);
                   end;
               end;
               admittance(node,node)=y_temp;
end;
for i=1:mn(1);
              admittance(put_in(i,1),put_in(i,2))=-put_in(i,3);
end;
admittance=conj(admittance')+admittance-tril(admittance,0);
voltage=linspace(1+0i,1+0i,sum_node);
voltage(1)=balance_power_node;
star_node=2;
end_node=sum_node;                      %由于没有涉及PV节点，暂时用此式表示
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------------------------牛顿拉夫逊算法实现------------------------------%
for mm=1:max_n;
%注入功率计算
e=real(voltage);
f=imag(voltage);
G=real(admittance);
B=imag(admittance);
for i=1:end_node;
        for j=1:end_node ;
            temp_P(i,j)=e(i)*(G(i,j)*e(j)-B(i,j)*f(j))+f(i)*(G(i,j)*f(j)+B(i,j)*e(j));
            temp_Q(i,j)=f(i)*(G(i,j)*e(j)-B(i,j)*f(j))-e(i)*(G(i,j)*f(j)+B(i,j)*e(j));
        end;
        P(i)=sum(temp_P(i,:));
        Q(i)=sum(temp_Q(i,:));
end;
%注入功率误差计算
power_n=P+Q*(1i);
d_power=power-power_n                %功率不平衡量
d_power_real=real(d_power);
d_power_imag=imag(d_power);
%注入电流计算
voltage_conj=conj(voltage);
power_n_conj=conj(power_n);
current=power_n_conj./voltage_conj;
%雅克比修正矩阵
a=real(current);
b=imag(current);
k=1;
for i=star_node:end_node;
               d_s_T(k:k+1,1)=[d_power_real(i);d_power_imag(i)];
               voltage_T(k:k+1,1)=[f(i);e(i)];                       k=k+2;
               m=(i-star_node)*2+1;
               for j=star_node:end_node;
                   if i==j
                       H(i,i)=-B(i,i)*e(i)+G(i,i)*f(i)+b(i);
                       N(i,i)= G(i,i)*e(i)+B(i,i)*f(i)+a(i);
                       J(i,i)=-G(i,i)*e(i)-B(i,i)*f(i)+a(i);
                       L(i,i)=-B(i,i)*e(i)+G(i,i)*f(i)-b(i);
                       R(i,i)= 2*f(i);S(i,i)= 2*e(i);
                   else
                       H(i,j)=-B(i,j)*e(i)+G(i,j)*f(i);
                       N(i,j)= G(i,j)*e(i)+B(i,j)*f(i);
                       J(i,j)=-N(i,j);R(i,j)= 0;
                       L(i,j)= H(i,j);S(i,j)= 0;           
                   end;
                   n=(j-star_node)*2+1;
                   Jacobian(m:(m+1),n:(n+1))=[H(i,j),N(i,j);J(i,j),L(i,j)];
               end;
end;
diag(Jacobian)';%雅克比对角量
%修正电压初值
J_1=inv(Jacobian);
d_v_T=J_1*d_s_T;
voltage_T=voltage_T+d_v_T;
k=1;
s=2;
temp=voltage;
for i=star_node:end_node
               voltage(s)=(voltage_T(k,1)*(1i)+voltage_T(k+1,1));
               k=k+2;
               s=s+1;
end
voltage ;%修正后的电压
voltage-temp;%电压修正量
if             max(d_v_T)<=error
               break;
end;
end;
%-----------------------------报告生成模块--------------------------------%
power_n(1)=P(1)+Q(1)*(1i);                        %平衡点的功率
for i=1:sum_node                                        %线路功率损耗
               for j=1:sum_node
power_ware(i,j)=-voltage(i)*conj((voltage(i)-voltage(j))*admittance(i,j));
               end;
end;
sum_power_n=sum(power_n);                         %总的损耗功率
temp=real(power_n);                                        %网络输电效率
efficiency=-sum(real(power_n(find(temp<0))))/sum(real(power_n(find(temp>0))))*100;
voltage_angle=angle(voltage)*180/pi;              %电压用指数形式表示
voltage_abs=abs(voltage);
voltage_e=[voltage_abs;voltage_angle];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%          最终结果             %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc                           %%%%%%%                
power_n                       %%%%%%%    日期：  2006.12.20
power_ware                    %%%%%%%    第一版
sum_power_n                   %%%%%%%          
efficiency                    %%%%%%%    http://hi.baidu.com/lyz_4534
voltage_e                     %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%