www.pudn.com > OPF.zip > Jacobi.m, change:2013-01-24,size:4278b


%% *************************************************************** 
%             作者:陆镛 
%             原创日期:2012年11月26号 
%             修改日期:2012年12月2号 
%             函数说明:修正方程中的不等式约束雅可比矩阵。 
%% *************************************************************** 
function [InequalityJmatrix,EqualityJmatrix,Jmatrix,Ycos,Ysin] = Jacobi(PvNum,... 
    PvI,GenNum,GenI,GenA,NodeNum,y,L,U,Z,W,AngleU,AmplitudeU,AmplitudeY,AngleY,I,J) 
%% 变量说明:****************************************************** 
% 输入参数说明:PvNum:Pv节点数  Pv:PV节点号  GenNum:发电机节点数  GenI:发电机节点号 
% GenA:发电机耗量特性二次项系数  NodeNum:节点数  y:拉格朗日乘子 Z、W:拉格朗日乘子对角阵 
% L、U:松弛变量对角阵  AngleU:电压相角   AmplitudeU:电压幅值  AmplitudeY:节点导纳矩阵幅值 
% AngleY:节点导纳矩阵相角 
% 输出参数说明:InequalityJmatrix:不等式约束矩阵  EqualityJmatrix:等式约束矩阵 
% Jmatrix:海森矩阵和不等式约束矩阵组成的雅可比矩阵 Ycos:导纳矩阵幅值点乘AngleUij减去导纳矩阵相角的余弦值 
% Ysin:导纳矩阵幅值点乘AngleUij减去导纳矩阵相角的正弦值 
%% 等式约束 
DiagU=sparse(1:NodeNum,1:NodeNum,AmplitudeU,NodeNum,NodeNum);             
AngleUij=sparse([I;J],[J;I],[AngleU(I)-AngleU(J);AngleU(J)-AngleU(I)],...    
    NodeNum,NodeNum);                                                     
Ysin=AmplitudeY.*spfun(@sin,AngleUij-AngleY);                           
Ycos=AmplitudeY.*spfun(@cos,AngleUij-AngleY);                               
                                                   
dhx=[DiagU*(Ysin*DiagU-sparse(1:NodeNum,1:NodeNum,Ysin*AmplitudeU,NodeNum,NodeNum)),... 
    sparse(1:NodeNum,1:NodeNum,Ycos*AmplitudeU,NodeNum,NodeNum)+DiagU*Ycos;... 
    -DiagU*(Ycos*DiagU-sparse(1:NodeNum,1:NodeNum,Ycos*AmplitudeU,NodeNum,NodeNum)),... 
    sparse(1:NodeNum,1:NodeNum,Ysin*AmplitudeU,NodeNum,NodeNum)+DiagU*Ysin];                                                         
EqualityJmatrix=[sparse(1:GenNum,GenI,-1,GenNum,2*NodeNum);sparse(1:PvNum... 
    ,PvI+NodeNum,-1,PvNum,2*NodeNum);dhx']; 
%% 不等式约束 
InequalityJmatrix=[sparse(1:GenNum+PvNum,1:GenNum+PvNum,1,GenNum+PvNum,... 
    PvNum+GenNum+NodeNum);[sparse(2*NodeNum,GenNum+PvNum) sparse(1+... 
    NodeNum:2*NodeNum,1:NodeNum,1,2*NodeNum,NodeNum)]]; 
%% 海森阵 
DiagU=sparse(1:NodeNum,1:NodeNum,AmplitudeU,NodeNum,NodeNum); 
yp=y(1:NodeNum); 
yq=y(1+NodeNum:2*NodeNum); 
DiagYp=sparse(1:NodeNum,1:NodeNum,yp,NodeNum,NodeNum); 
DiagYq=sparse(1:NodeNum,1:NodeNum,yq,NodeNum,NodeNum); 
 
d2hxy3H1=(-sparse(1:NodeNum,1:NodeNum,Ycos'*DiagU*yp,NodeNum,NodeNum)+DiagU*... 
    DiagYp*Ycos)*DiagU-(sparse(1:NodeNum,1:NodeNum,Ycos*AmplitudeU)-... 
    DiagU*Ycos')*DiagU*DiagYp; 
d2hxy3H2=(-sparse(1:NodeNum,1:NodeNum,Ysin*AmplitudeU,NodeNum,NodeNum)+... 
    DiagU*Ysin')*DiagU*DiagYq-(-DiagU*DiagYq*Ysin+sparse(1:NodeNum,... 
    1:NodeNum,Ysin'*DiagU*yq,NodeNum,NodeNum))*DiagU; 
d2hxy3H=d2hxy3H1+d2hxy3H2; 
d2hxy3N1=(-sparse(1:NodeNum,1:NodeNum,Ysin*AmplitudeU,NodeNum,NodeNum)+... 
    DiagU*Ysin')*DiagYp+(-DiagU*DiagYp*Ysin+sparse(1:NodeNum,1:NodeNum,... 
    Ysin'*DiagU*yp,NodeNum,NodeNum)); 
d2hxy3N2=(sparse(1:NodeNum,1:NodeNum,Ycos*AmplitudeU,NodeNum,NodeNum)-... 
    DiagU*Ycos')*DiagYq+DiagU*DiagYq*Ycos-sparse(1:NodeNum,1:NodeNum,... 
    Ycos'*DiagU*yq,NodeNum,NodeNum); 
d2hxy3N=d2hxy3N1+d2hxy3N2; 
d2hxy3J1=sparse(1:NodeNum,1:NodeNum,Ysin'*DiagU*yp,NodeNum,NodeNum)+... 
    DiagYp*Ysin*DiagU-DiagYp*sparse(1:NodeNum,1:NodeNum,Ysin*AmplitudeU... 
    ,NodeNum,NodeNum)-Ysin'*DiagU*DiagYp; 
d2hxy3J2=Ycos'*DiagU*DiagYq+DiagYq*sparse(1:NodeNum,1:NodeNum,... 
    Ycos*AmplitudeU,NodeNum,NodeNum)-sparse(1:NodeNum,1:NodeNum,... 
    Ycos'*DiagU*yq,NodeNum,NodeNum)-DiagYq*Ycos*DiagU; 
d2hxy3J=d2hxy3J1+d2hxy3J2; 
d2hxy3L1=Ycos'*DiagYp+DiagYp*Ycos; 
d2hxy3L2=Ysin'*DiagYq+DiagYq*Ysin; 
d2hxy3L=d2hxy3L1+d2hxy3L2; 
d2hxy3=[d2hxy3H d2hxy3N;d2hxy3J d2hxy3L]; 
 
d2hxy=[sparse(PvNum+GenNum,PvNum+GenNum+2*NodeNum);sparse(2*NodeNum... 
    ,PvNum+GenNum),d2hxy3]; 
%%  
HessenMatrix=-sparse(1:GenNum,1:GenNum,2*GenA,PvNum+GenNum+2*NodeNum,PvNum+GenNum... 
    +2*NodeNum)+d2hxy-InequalityJmatrix*(L\Z-U\W)*InequalityJmatrix'; 
%% 由海森阵和不等式约束矩阵组成的雅可比矩阵 
Jmatrix=[HessenMatrix,EqualityJmatrix;EqualityJmatrix',sparse(2*NodeNum,2*NodeNum)]; 
end