最优潮流Matlab程序

function [ ] = optflow（  ）
%UNtitleD Summary of this function goes here
%   Detailed explanation goes here
y=[2.1569-8.6274i          0         -0.5882+2.3529i    -0.3922+1.5686i   -1.1765+4.7059i;
         0          2.3530-9.4118i   -1.1765+4.7059i           0          -1.1765+4.7059i;
  -0.5882+2.3529i   -1.1765+4.7059i   2.3529-9.4117i     -0.5882+2.3529i         0;
  -0.3922+1.5686i          0         -0.5882+2.3529i     0.9804-3.9215i          0;
  -1.1765+4.7059i   -1.1765+4.7059i         0                  0           2.3530-9.4118i]  ;%导纳矩阵
g=real（y）;
b=imag（y）;
u（:1）=[1;1;1;1;1];        %各节点初始电压值
ct（:1）=[0;0;0;0;0];       %初始角度
l=ones（114）;
uu=ones（114）;
z=ones（114）;
w=-0.5*ones（114）;
y=ones（110）;
y（1:2:9）=（1e-10）*y（1:2:9）;
y（2:2:10）=（-1e-10）*y（2:2:10）;
p0=[1;0;-0.6;-0.4;-0.6];   %各节点的输入有功功率
q0=[0.5;0;-0.3;-0.1;-0.2]; %各节点的输入无功功率
hPG=zeros（210）;
hPG（14*2-1）=-1;
hPG（25*2-1）=-1;%发电机1在4节点上，发电机2在5节点上
hQR（14*2）=-1;
hQR（25*2）=-1;
hx=zeros（1010）;
hxx=zeros（1410）;
gxx=zeros（1414）;
k=1;
 for i=1:5
    for j=1:5
        if i~=j
            hh（ij）=-u（ik）*u（jk）*（g（ij）*sin（ct（ik）-ct（jk））-b（ij）*cos（ct（ik）-ct（jk）））;
            jj（ij）=-u（ik）*u（jk）*（g（ij）*cos（ct（ik）-ct（jk））+b（ij）*sin（ct（ik）-ct（jk）））;
            nn（ij）=-u（ik）*（g（ij）*cos（ct（ik）-ct（jk））+b（ij）*sin（ct（ik）-ct（jk）））;