基于粒子群算法的无功优化MATLAB源代码，IEEE30节点

function [SumQg,Only_PLoss] = IPM_MainCalculation(lanbda) % 控制变量只有发电机机端电压 % 改动GradsMatrix.m，GradsHeisenMatrix.m，PLoss_Calculation.m clc ND_COUNT = 50;% 内点算法迭代次数 CENTER_PARAM = 0.1;% 向心参数(电力系统中一般为0.05-0.2) STEP_UNIT = 0.9995;% 补偿单位（安全因子，保证某些数值为正） DEF_PRESICION = 0.000001;%收敛精度 %% ************************************** 数据初始化 **************************************** [TransFormer_Branch, Normal_Branch, PQ_Node, PV_Node, Swing_Node, Node_Num] = RE_IEEE14_data; % 输入节点数据，Swing_Node为平衡节点 [ m_Node, Total_PLoad,Total_QLoad, Node_Num, GenCount, nVaribles, nEquals, nNEquals ] = Form_IPM_data14( PQ_Node, PV_Node, Swing_Node, Node_Num ); [GB, m_G, m_B] = FormGB_former(TransFormer_Branch, Normal_Branch, PQ_Node,PV_Node, Node_Num, m_Node); % GB为导纳矩阵，m_G为其实部，m_G为其虚部，变压器支路已处理 [GB, G, B] = FormGB_former(TransFormer_Branch, Normal_Branch, PQ_Node, PV_Node, Node_Num, m_Node);% [ PLoss,voltagee,angle ] = Power_Calculation( G, B, PQ_Node, PV_Node, Swing_Node, Node_Num, m_Node); PLoss; [PLoss,SumQg,Only_PLoss] = PLoss_Calculation( m_Node, Node_Num, m_G, m_B, lanbda ) % PQ节点 Type = 0 % PV节点 Type = 1 % 平衡节点 Type = 2 % 两次算得的Ploss不同，因为节点电压不同 %% ************************************* 程序初始化 **************************************** m_l = zeros ( nNEquals, 1 ); m_u = zeros ( nNEquals, 1 ); % 松弛变量 j = 1; for i = 1:Node_Num if ( m_Node(i).Type == 1 ) m_l (j) = m_Node(i).Qg - m_Node(i).QDnLmt; m_u (j) = - m_Node(i).Qg + m_Node(i).QUpLmt; a(j) = m_l(j); b(j) = m_u(j); j = j + 1; end end for i = 1:Node_Num m_l (j) = m_Node(i).Vol - m_Node(i).VolDnLmt; m_u (j) = - m_Node(i).Vol + m_Node(i).VolUpLmt; a(j) = m_l(j); b(j) = m_u(j); j = j + 1; end m_z = ones( nNEquals, 1 ); m_w = -0.5 * ones( nNEquals, 1 ); m_y = -0.01 * ones( nEquals, 1 ); %% ****************************迭代程序计算****************************** flag = false; nCount = 0; while ( ( ~flag ) & ( nCount <= ND_COUNT ) ) nCount = nCount + 1; % 计算补偿间隙 fGap = m_l'*m_z - m_u'*m_w; Gap( nCount ) = fGap; if ( fGap < DEF_PRESICION ) flag = true; break; end % 计算障碍因子 m_fR = CENTER_PARAM * fGap / 2 / nNEquals; % 修正各节点的值 [ m_Node ] = ModifyNodeValue( m_Node, m_G, m_B, Node_Num ); % 计算等约束的雅克比矩阵 Jocabi_H = zeros ( nVaribles, nEquals ); Jocabi_H = EqualJacobiMatrix( m_Node, m_G, m_B, Node_Num, GenCount ); % 计算等约束的海森矩阵 Hession_H = zeros ( nVaribles ); Hession_H = EqualHeisenMatrix( m_Node, m_G, m_B, m_y, GenCount, Node_Num, nEquals, nVaribles ); % 计算不等约束的雅克比矩阵 Jocabi_G = zeros ( nVaribles, nNEquals ); Jocabi_G = InequalJacobiMatrix( GenCount, Node_Num); %% 计算目标函数的梯度向量,在这部分修改 Grads_obj = zeros ( nVaribles, 1 ); Grads_obj = GradsMatrix( m_Node, m_G, m_B, GenCount, Node_Num,lanbda); % 计算目标函数的海森矩阵,在这部分修改 Hession_obj = zeros ( nVaribles ); Hession_obj = GradsHeisenMatrix( m_Node, m_G, m_B, GenCount, Node_Num,lanbda); %% 计算△x和△y L = zeros ( nNEquals ); Z = zeros ( nNEquals ); U = zeros ( nNEquals ); W = zeros ( nNEquals ); for i = 1:nNEquals L( i, i ) = m_l(i); U( i, i ) = m_u(i); Z( i, i ) = m_z(i); W( i, i ) = m_w(i); end %% ***************** 求解一阶最优性条件解 ********************************* L_x = Grads_obj - Jocabi_H * m_y - Jocabi_G * ( m_z + m_w ); L_y = zeros ( nEquals, 1); k = 1; for i = 1:Node_Num if ( m_Node(i).Type ~= 2 ) L_y(k) = m_Node(i).Pg - m_Node(i).Pl - m_Node(i).P; k = k+1; end end for i = 1:Node_Num L_y(k) = m_Node(i).Qg - m_Node(i).Ql - m_Node(i).Q; k = k+1; end Power_err( nCount ) = abs( max (L_y) ); L_z = zeros ( nNEquals, 1 ); L_w = zeros ( nNEquals, 1 ); k = 1; for i = 1:Node_Num if ( m_Node(i).Type == 1 ) L_z(k) = m_Node(i).Qg - m_Node(i).QDnLmt - m_l(k); L_w(k) = m_Node(i).Qg - m_Node(i).QUpLmt + m_u(k); k = k+1; end end for i = 1:Node_Num L_z(k) = m_Node(i).Vol - m_Node(i).VolDnLmt - m_l(k); L_w(k) = m_Node(i).Vol - m_Node(i).VolUpLmt + m_u(k); k = k+1; end L_l = L * Z * ones( nNEquals, 1 ) - m_fR * ones( nNEquals, 1 ); L_u = U * W * ones( nNEquals, 1 ) + m_fR * ones( nNEquals, 1 ); %% ***************** 求解Δx、Δy、Δz、Δl、Δw、Δu变量 *************************************** Hession = Hession_H - Hession_obj; %不等式的海森矩阵为零 Hession = Hession - Jocabi_G * ( inv(L)*Z - inv(U)*W ) * Jocabi_G.'; %disp ( Hession ) HH = [ Hession, Jocabi_H Jocabi_H.', zeros(nEquals) ]; L_xx = L_x + Jocabi_G * ( inv(L)*( L_l + Z * L_z ) + inv(U) * ( L_u - W * L_w) ); xy_right = [ L_xx; - L_y]; xy =( HH ) \ xy_right; delt_x = xy ( 1: nVaribles, 1); delt_y = xy ( nVaribles+1 : nVaribles+nEquals, 1 ); delt_u = -L_w - Jocabi_G.' * delt_x; delt_w = - inv( U ) * ( L_u + W * delt_u); delt_l = L_z - L_w - delt_u; delt_z = -inv(L) *( L_l + Z * delt_l ); %% ************************** 修正各变量值 ****************************************** % 计算原对偶步长：Tp, Td fp = 1.0; fd = 1.0; for i = 1: nNEquals f = - m_l(i) / delt_l(i); f1(i) = - m_l(i) / delt_l(i); if ( f > 0 && f < fp ) fp = f; end end f1; fp; for i = 1: nNEquals f = - m_u(i) / delt_u(i); f2(i) = - m_u(i) / delt_u(i); if ( f > 0 && f < fp ) fp = f; end end f2; fp; for i = 1:nNEquals f = - m_z(i) / delt_z(i); f3(i) = - m_z(i) / delt_z(i); if ( f > 0 && f < fd ) fd = f; end end f3; fd; for i = 1:nNEquals f = - m_w(i) / delt_w(i); f4(i) = - m_w(i) / delt_w(i); if ( f > 0 && f < fd ) fd = f; end end f4; fd; fp = fp * STEP_UNIT; fd = fd * STEP_UNIT; %修正x k = 1; for i = 1:Node_Num if ( m_Node(i).Type == 1 ) m_Node(i).Qg = m_Node(i).Qg + fp * delt_x(k); Qg( k ) = m_Node(i).Qg; k = k + 1; end end for i = 1:Node_Num m_Node(i).Vol = m_Node(i).Vol + fp * delt_x(k); vol( i ) = m_Node(i).Vol; k = k + 1; end for i = 1:Node_Num if ( m_Node(i).Type ~= 2 ) m_Node(i).Angle = m_Node(i).Angle + fp * delt_x(k); k = k + 1; end angle(i) = m_Node(i).Angle; end m_l = m_l + fp * delt_l; m_u = m_u + fp * delt_u; m_y = m_y + fd * delt_y; m_z = m_z + fd * delt_z; m_w = m_w + fd * delt_w; [PLoss,SumQg,Only_PLoss] = PLoss_Calculation( m_Node, Node_Num, m_G, m_B, lanbda); P_PLoss ( nCount ) = PLoss; vol; angle*180/3.14; for i = 1:Node_Num q1(i) = m_Node(i).VolUpLmt - vol(i); w1(i) = vol(i) - m_Node(i).VolDnLmt; end q1; w1; end %% %nCount = nCount %flag P_slack = 0; if ( flag == true ) %vol %angle*180/3.14 %P_PLoss Qg %Total_PLoad %Total_QLoad for i = 1:Node_Num f = m_Node(Node_Num).Angle - m_Node(i).Angle; P_slack = P_slack + m_Node(Node_Num).Vol * m_Node(i).Vol * ( m_G(Node_Num,i)*cos(f) + m_B(Node_Num,i)*sin(f) ); end P_slack = P_slack + m_Node( Node_Num ).Pl; QLoss = 0.0; for i = 1:Node_Num for j = 1:Node_Num f = m_Node(i).