facat
/
stateestimation-ipm
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

似乎公式没推错，只是牛顿法对初值敏感，带入潮流的电压和相角迭代几次就可以收敛。 

Signed-off-by: dmy <dugg@21cn.com>
This commit is contained in:
dmy 2013-09-26 22:04:35 +08:00
parent a9d49d8edd
commit 59f7810cab
14 changed files with 484 additions and 68 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
.gitignore vendored
View File

@ -1,6 +1,7 @@
/Powerflow
/Powerflow
/公式/*.aux
/公式/*.log
/公式/*.pdf
/公式/*.synctex.gz(busy)
*.asv
*.asv
*.mat

2
BranchI.m
View File

@ -1,5 +1,5 @@
function [ output_args ] = BranchI( V,lineI,lineJ,lineR,lineX,lineB2 )
output_args=(V(lineI)-V(lineJ))./(lineR+1j*lineX)+V(lineI).*lineB2;
output_args=(V(lineI)-V(lineJ))./(lineR+1j*lineX)+V(lineI).*(1j*lineB2);
end

1
Powerflow Submodule

@ -0,0 +1 @@
Subproject commit a8da7f966d23827f01fec553a52bab3186c8a35c

BIN
mBranchI.mat
View File

Binary file not shown.

BIN
mBranchP.mat
View File

Binary file not shown.

BIN
mBranchQ.mat
View File

Binary file not shown.

BIN
mPG.mat
View File

Binary file not shown.

BIN
mQG.mat
View File

Binary file not shown.

BIN
mTransP.mat
View File

Binary file not shown.

BIN
mTransQ.mat
View File

Binary file not shown.

BIN
mVolt.mat
View File

Binary file not shown.

461
run.m
View File

@ -2,7 +2,11 @@ clear
clc
% yalmip('clear')
addpath('.\Powerflow')
[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee14.dat', '0');
[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('s1047.dat', '0');
%% 节点总数
busNum=length(Volt);
%% 变量总数
varNum=busNum*2;
%% 开始生成量测量
sigma=0.03;% 标准差
%% 电压
@ -33,12 +37,16 @@ mBranchI=rBranchI.*(normrnd(0,sigma,length(rBranchI),1)+1);%֧·
%% 支路功率
rBranchP= BranchP( cmpV,cmpBranchI,lineI,lineB2 );
mBranchP=rBranchP.*(normrnd(0,sigma,length(rBranchP),1)+1);%支路功率量测量
%load('mBranchP');
%mBranchP=[0.284828008908987,-0.0282291720369099,-0.296295663000711]';
rBranchQ=BranchQ( cmpV,cmpBranchI,lineI,lineB2 );
mBranchQ=rBranchQ.*(normrnd(0,sigma,length(rBranchQ),1)+1);%支路功率量测量
%% 变压器功率
transI=newwordParameter.trans.transI;
transJ=newwordParameter.trans.transJ;
% newwordParameter.trans.transK=ones(length(newwordParameter.trans.transK),1);
transK=newwordParameter.trans.transK;
% transK=ones(length(transK),1);
transR=newwordParameter.trans.transR;
transX=newwordParameter.trans.transX;
transG=real(1./(transR+1j*transX));
@ -73,7 +81,7 @@ onlyQG=setdiff(QGi,PDQDi);
measureSigma=abs(([rVolt;rBranchP;rBranchQ;rTransP;rTransQ].*sigma));
measureSigma(measureSigma<1e-6)=mean(measureSigma(measureSigma>1e-6));
W=sparse(diag(1./measureSigma.^2)) ;
% W=sparse(1:length(W),1:length(W),400,length(W),length(W));
W=sparse(1:length(W),1:length(W),1,length(W),length(W));
%% 冗余度计算
stateVarCount=2*length(Volt);
measurements=length(mVolt)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG)+length(mTransP)+length(mTransQ);
@ -99,8 +107,17 @@ fprintf('
%% 自己写的微分代码
% 初始化一些值
SEVolt=sparse(ones(length(mVolt),1));
SEVolt(Balance)=rVolt(Balance);
SEVAngle=sparse(zeros(length(mVolt),1));
SEVolt=rVolt;
% SEVolt(1:2:end)=1.2;
% load('SEVolt');
% load('SEVAngle');
%SEVolt(Balance)=rVolt(Balance);
SEVAngle=-0.1*sparse(ones(length(mVolt),1));
SEVAngle=rVAngel;
% SEVolt2= load('SEVolt2');
% SEVolt2=SEVolt2.SEVolt;
% SEVAngle2= load('SEVAngle2');
% SEVAngle2=SEVAngle2.SEVAngle;
maxD=1000;
Iteration=0;
optimalCondition=100;
@ -110,8 +127,8 @@ v=2;
ojbFunDecrease=1000;% 目标函数下降
% 以下都是Jacobi矩阵
% while max(abs(g))>1e-5;
% while maxD>1e-5
while max(abs(optimalCondition))>eps
while maxD>1e-5
%while max(abs(optimalCondition))>eps
% 电压
dV_dV=sparse(1:length(mVolt),1:length(mVolt),1,length(mVolt),length(mVolt));%电压量测量的微分
dV_dTyta=sparse(length(mVolt),length(mVolt));
@ -138,6 +155,7 @@ while max(abs(optimalCondition))>eps
)...
,length(lineI),length(mVolt));%线路的
dLPij=[dLPij_dVi+dLPij_dVj,dLPij_dThetai+dLPij_dThetaj];
dLQij_dVi=sparse(1:length(lineI),lineI, ...
-SEVolt(lineJ).*( ...
lineG.*sin(SEVAngle(lineI)-SEVAngle(lineJ)) -lineB.*cos(SEVAngle(lineI)-SEVAngle(lineJ))...
@ -159,6 +177,7 @@ while max(abs(optimalCondition))>eps
lineG.*cos(SEVAngle(lineI)-SEVAngle(lineJ)) +lineB.*sin(SEVAngle(lineI)-SEVAngle(lineJ))...
) ...
,length(lineI),length(mVolt));%线路的
dLQij=[dLQij_dVi+dLQij_dVj,dLQij_dThetai+dLQij_dThetaj];
% 变压器支路
dTPij_dVi=sparse(1:length(transI),transI, ...
-SEVolt(transJ)./transK.*( ...
@ -232,84 +251,402 @@ while max(abs(optimalCondition))>eps
% nodePQ=[nodeP;nodeQ];
% c=nodePQ(zerosInjectionIndex);
%% 进入迭代
% 一阶导数
H=[dV_dV,dV_dTyta;
dLPij_dVi+dLPij_dVj,dLPij_dThetai+dLPij_dThetaj ;
dLQij_dVi+dLQij_dVj,dLQij_dThetai+dLQij_dThetaj ;
dTPij_dVi+dTPij_dVj,dTPij_dThetai+dTPij_dThetaj;
dTQij_dVi+dTQij_dVj,dTQij_dThetai+dTQij_dThetaj];%jacobi
% 二阶导数
SEBranchI=BranchI( SEVolt.*exp(1j*SEVAngle),lineI,lineJ,lineR,lineX,lineB2 );%复数支路电流
SEBranchP=BranchP( SEVolt.*exp(1j*SEVAngle),SEBranchI,lineI,lineB2 );
SEBranchQ=BranchQ( SEVolt.*exp(1j*SEVAngle),SEBranchI,lineI,lineB2 );
SETransP=TransPower( newwordParameter,SEVolt,SEVAngle );
SETransQ=TransReactivePower( newwordParameter,SEVolt,SEVAngle );
h=[SEVolt;SEBranchP;SEBranchQ;SETransP;SETransQ];
% load('mBranchP');
% mBranchQ=mBranchP;
% mTransP=mBranchP(1:length(mTransP));
% mTransQ=mBranchP(1:length(mTransQ));
% mVolt=mBranchP(1:length(mVolt));
%mBranchP=[0.284828008908987,-0.0282291720369099,-0.296295663000711]';
% load('mTransP');
% mTransP=mTransP';
% mTransQ=mTransP;
z=[mVolt;mBranchP;mBranchQ;mTransP;mTransQ];
G=H'*W*H;
g=-H'*W*(z-h);
Whz=W*(h-z);%形成海森矩阵时使用
% 形成大的求解矩阵
% a=[G C';
% C zeros(size(C,1),size(C,1))];
% b=[-g;-c];
% 利用 Levenber-Marquardt
% if Iteration==0
% mu=max(diag(G));
% end
%% 一阶导数
H=[dV_dV,dV_dTyta;
dLPij_dVi+dLPij_dVj,dLPij_dThetai+dLPij_dThetaj ;
dLQij_dVi+dLQij_dVj,dLQij_dThetai+dLQij_dThetaj ;
dTPij_dVi+dTPij_dVj,dTPij_dThetai+dTPij_dThetaj;
dTQij_dVi+dTQij_dVj,dTQij_dThetai+dTQij_dThetaj];%jacobi
%% 二阶导数
% 线路二阶导数
% 对电压求导
d2LPij_dVi2=sparse(lineI,lineI, ...
2*lineG.*Whz(length(mVolt)+1:length(mVolt)+length(lineI)),...
varNum,varNum);%线路有功海森矩阵
d2LPij_dVidVj=sparse(lineI,lineJ, ...
Whz(length(mVolt)+1:length(mVolt)+length(lineI))...
.* ( -(cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG + sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB)), ...
varNum,varNum);%线路有功海森矩阵
d2LPij_dVj2=0;%线路有功海森矩阵
%%
% d2LPij_dVi2=sparse(lineI,lineI, ...
% 2*lineG ,...
% varNum,varNum);%线路有功海森矩阵
% d2LPij_dVidVj=sparse(lineI,lineJ, ...
% 1 ...
% .* ( -(cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG + sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB)), ...
% varNum,varNum);%线路有功海森矩阵
% d2LPij_dVj2=0;%线路有功海森矩阵
% t1=d2LPij_dVi2+d2LPij_dVidVj+d2LPij_dVidVj';
% 对相角求导
d2LPij_dTi2=sparse(lineI+busNum,lineI+busNum, ...
Whz(length(mVolt)+1:length(mVolt)+length(lineI))...
.*SEVolt(lineI).*SEVolt(lineJ).*( cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG + sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB), ...
varNum,varNum);%线路有功海森矩阵
d2LPij_dTidTj=sparse(lineI+busNum,lineJ+busNum, ...
Whz(length(mVolt)+1:length(mVolt)+length(lineI))...
.*(-SEVolt(lineI).*SEVolt(lineJ).*( cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG + sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB) ), ...
varNum,varNum);%线路有功海森矩阵
d2LPij_dTj2=sparse(lineJ+busNum,lineJ+busNum, ...
Whz(length(mVolt)+1:length(mVolt)+length(lineI))...
.*SEVolt(lineI).*SEVolt(lineJ).*( cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG + sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB), ...
varNum,varNum);%线路有功海森矩阵
%%
% d2LPij_dTi2=sparse(lineI+busNum,lineI+busNum, ...
% 1 ...
% .*SEVolt(lineI).*SEVolt(lineJ).*( cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG + sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB), ...
% varNum,varNum);%线路有功海森矩阵
% d2LPij_dTidTj=sparse(lineI+busNum,lineJ+busNum, ...
% 1 ...
% .*(-SEVolt(lineI).*SEVolt(lineJ).*( cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG + sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB) ), ...
% varNum,varNum);%线路有功海森矩阵
% d2LPij_dTj2=sparse(lineJ+busNum,lineJ+busNum, ...
% 1 ...
% .*SEVolt(lineI).*SEVolt(lineJ).*( cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG + sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB), ...
% varNum,varNum);%线路有功海森矩阵
% t2=d2LPij_dTi2+d2LPij_dTidTj+d2LPij_dTj2+d2LPij_dTidTj';
%%
% d2LPij_dTi2=sparse(lineI+busNum,lineI+busNum, ...
% 1 ...
% .*SEVolt(lineI).*SEVolt(lineJ).*( cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG + sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB), ...
% varNum,varNum);%线路有功海森矩阵
% d2LPij_dTidTj=sparse(lineI+busNum,lineJ+busNum, ...
% 1 ...
% .*(-SEVolt(lineI).*SEVolt(lineJ).*( cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG + sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB) ), ...
% varNum,varNum);%线路有功海森矩阵
% d2LPij_dTj2=sparse(lineJ+busNum,lineJ+busNum, ...
% 1 ...
% .*SEVolt(lineI).*SEVolt(lineJ).*( cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG + sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB), ...
% varNum,varNum);%线路有功海森矩阵
% t=d2LPij_dTi2+d2LPij_dTidTj+d2LPij_dTidTj'+d2LPij_dTj2;
%%
% 对电压相角求导
d2LPij_dVidTi=sparse(lineI,lineI+busNum, ...
Whz(length(mVolt)+1:length(mVolt)+length(lineI))...
.* (SEVolt(lineJ).* (sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG - cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
varNum,varNum);%线路有功海森矩阵
d2LPij_dVidTj=sparse(lineI,lineJ+busNum, ...
Whz(length(mVolt)+1:length(mVolt)+length(lineI))...
.* (-SEVolt(lineJ).* (sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG - cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
varNum,varNum);%线路有功海森矩阵
d2LPij_dVjdTi=sparse(lineJ,lineI+busNum, ...
Whz(length(mVolt)+1:length(mVolt)+length(lineI))...
.* (SEVolt(lineI).* (sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG - cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
varNum,varNum);%线路有功海森矩阵
d2LPij_dVjdTj=sparse(lineJ,lineJ+busNum, ...
Whz(length(mVolt)+1:length(mVolt)+length(lineI))...
.* (-SEVolt(lineI).* (sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG - cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
varNum,varNum);%线路有功海森矩阵
%%
% d2LPij_dVidTi=sparse(lineI,lineI+busNum, ...
% 1 ...
% .* (SEVolt(lineJ).* (sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG - cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
% varNum,varNum);%线路有功海森矩阵
% d2LPij_dVidTj=sparse(lineI,lineJ+busNum, ...
% 1 ...
% .* (-SEVolt(lineJ).* (sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG - cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
% varNum,varNum);%线路有功海森矩阵
% d2LPij_dVjdTi=sparse(lineJ,lineI+busNum, ...
% 1 ...
% .* (SEVolt(lineI).* (sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG - cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
% varNum,varNum);%线路有功海森矩阵
% d2LPij_dVjdTj=sparse(lineJ,lineJ+busNum, ...
% 1 ...
% .* (-SEVolt(lineI).* (sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG - cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
% varNum,varNum);%线路有功海森矩阵
% t3=d2LPij_dVidTi+d2LPij_dVidTj+d2LPij_dVjdTi+d2LPij_dVjdTj;
% t=t1+t2+t3+t3';
% 对电压求导
d2LQij_dVi2=sparse(lineI,lineI, ...
-2*(lineB+lineB2).*Whz(length(mVolt)+length(lineI)+1:length(mVolt)+length(lineI)*2),...
varNum,varNum);%线路无功海森矩阵
d2LQij_dVidVj=sparse(lineI,lineJ, ...
Whz(length(mVolt)+length(lineI)+1:length(mVolt)+length(lineI)*2) ...
.*( -(sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG - cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB)), ...
varNum,varNum);%线路无功海森矩阵
d2LQij_dVj2=0;%线路无功海森矩阵
% d2LQij_dVi2=sparse(lineI,lineI, ...
% -2*(lineB+lineB2),...
% varNum,varNum);%线路无功海森矩阵
% d2LQij_dVidVj=sparse(lineI,lineJ, ...
% 1 ...
% .*( -(sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG - cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB)), ...
% varNum,varNum);%线路无功海森矩阵
% d2LQij_dVj2=0;%线路无功海森矩阵
% t1=d2LQij_dVi2+d2LQij_dVidVj+d2LQij_dVj2+d2LQij_dVidVj';
% 对相角求导
d2LQij_dTi2=sparse(lineI+busNum,lineI+busNum, ...
Whz(length(mVolt)+length(lineI)+1:length(mVolt)+length(lineI)*2)...
.*SEVolt(lineI).*SEVolt(lineJ).*( sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG - cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB), ...
varNum,varNum);%线路无功海森矩阵
d2LQij_dTidTj=sparse(lineI+busNum,lineJ+busNum, ...
Whz(length(mVolt)+length(lineI)+1:length(mVolt)+length(lineI)*2)...
.*(-SEVolt(lineI).*SEVolt(lineJ).*( sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG - cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB)), ...
varNum,varNum);%线路无功海森矩阵
d2LQij_dTj2=sparse(lineJ+busNum,lineJ+busNum, ...
Whz(length(mVolt)+length(lineI)+1:length(mVolt)+length(lineI)*2)...
.*SEVolt(lineI).*SEVolt(lineJ).*( sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG - cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB), ...
varNum,varNum);%线路无功海森矩阵
% d2LQij_dTi2=sparse(lineI+busNum,lineI+busNum, ...
% 1 ...
% .*SEVolt(lineI).*SEVolt(lineJ).*( sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG - cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB), ...
% varNum,varNum);%线路无功海森矩阵
% d2LQij_dTidTj=sparse(lineI+busNum,lineJ+busNum, ...
% 1 ...
% .*(-SEVolt(lineI).*SEVolt(lineJ).*( sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG - cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB)), ...
% varNum,varNum);%线路无功海森矩阵
% d2LQij_dTj2=sparse(lineJ+busNum,lineJ+busNum, ...
% 1 ...
% .*SEVolt(lineI).*SEVolt(lineJ).*( sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG - cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB), ...
% varNum,varNum);%线路无功海森矩阵
% t2=d2LQij_dTi2+d2LQij_dTidTj+d2LQij_dTj2+d2LQij_dTidTj';
% 对电压相角求导
d2LQij_dVidTi=sparse(lineI,lineI+busNum, ...
Whz(length(mVolt)+length(lineI)+1:length(mVolt)+length(lineI)*2)...
.* (-SEVolt(lineJ).* (cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG + sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
varNum,varNum);%线路有功海森矩阵
d2LQij_dVidTj=sparse(lineI,lineJ+busNum, ...
Whz(length(mVolt)+length(lineI)+1:length(mVolt)+length(lineI)*2)...
.* (SEVolt(lineJ).* (cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG + sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
varNum,varNum);%线路有功海森矩阵
d2LQij_dVjdTi=sparse(lineJ,lineI+busNum, ...
Whz(length(mVolt)+length(lineI)+1:length(mVolt)+length(lineI)*2)...
.* (-SEVolt(lineI).* (cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG + sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
varNum,varNum);%线路有功海森矩阵
d2LQij_dVjdTj=sparse(lineJ,lineJ+busNum, ...
Whz(length(mVolt)+length(lineI)+1:length(mVolt)+length(lineI)*2)...
.* (SEVolt(lineI).* (cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG + sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
varNum,varNum);%线路有功海森矩阵
% d2LQij_dVidTi=sparse(lineI,lineI+busNum, ...
% 1 ...
% .* (-SEVolt(lineJ).* (cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG + sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
% varNum,varNum);%线路有功海森矩阵
% d2LQij_dVidTj=sparse(lineI,lineJ+busNum, ...
% 1 ...
% .* (SEVolt(lineJ).* (cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG + sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
% varNum,varNum);%线路有功海森矩阵
% d2LQij_dVjdTi=sparse(lineJ,lineI+busNum, ...
% 1 ...
% .* (-SEVolt(lineI).* (cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG + sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
% varNum,varNum);%线路有功海森矩阵
% d2LQij_dVjdTj=sparse(lineJ,lineJ+busNum, ...
% 1 ...
% .* (SEVolt(lineI).* (cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG + sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
% varNum,varNum);%线路有功海森矩阵
% t3=d2LQij_dVidTi+d2LQij_dVidTj+d2LQij_dVjdTi+d2LQij_dVjdTj;
% t=t1+t2+t3;
%% 变压器二阶导数
% 对电压求导
d2TPij_dVi2=sparse(transI,transI, ...
2*transG./transK./transK.*Whz(length(mVolt)+2*length(lineI)+1:length(mVolt)+2*length(lineI)+length(transI)),...
varNum,varNum);%变压器有功海森矩阵
d2TPij_dVidVj=sparse(transI,transJ, ...
Whz(length(mVolt)+2*length(lineI)+1:length(mVolt)+2*length(lineI)+length(transI))...
.* ( -(cos(SEVAngle(transI) -SEVAngle(transJ)).*transG + sin(SEVAngle(transI) -SEVAngle(transJ)).*transB)./transK), ...
varNum,varNum);%变压器有功海森矩阵
d2TPij_dVj2=0;%变压器有功海森矩阵
% d2TPij_dVi2=sparse(transI,transI, ...
% 2*transG./transK./transK.*1,...
% varNum,varNum);%变压器有功海森矩阵
% d2TPij_dVidVj=sparse(transI,transJ, ...
% 1 ...
% .* ( -(cos(SEVAngle(transI) -SEVAngle(transJ)).*transG + sin(SEVAngle(transI) -SEVAngle(transJ)).*transB)./transK), ...
% varNum,varNum);%变压器有功海森矩阵
% d2TPij_dVj2=0;%变压器有功海森矩阵
ggg=H(busNum+length(lineI)*2+1:busNum+length(lineI)*2+length(transI),:)'*H(busNum+length(lineI)*2+1:busNum+length(lineI)*2+length(transI),:)+d2TPij_dVi2+d2TPij_dVidVj+d2TPij_dVidVj';
% d2TPij_dVi2=sparse(transI,transI, ...
% 2*transG./transK./transK.*1,...
% varNum,varNum);%变压器有功海森矩阵
% d2TPij_dVidVj=sparse(transI,transJ, ...
% 1 ...
% .* ( -(cos(SEVAngle(transI) -SEVAngle(transJ)).*transG + sin(SEVAngle(transI) -SEVAngle(transJ)).*transB)./transK), ...
% varNum,varNum);%变压器有功海森矩阵
% d2TPij_dVj2=0;%变压器有功海森矩阵
% t1=d2TPij_dVi2+d2TPij_dVidVj+d2TPij_dVidVj';
% 对相角求导
d2TPij_dTi2=sparse(transI+busNum,transI+busNum, ...
Whz(length(mVolt)+2*length(lineI)+1:length(mVolt)+2*length(lineI)+length(transI))...
.*SEVolt(transI).*SEVolt(transJ)./transK.*( cos( SEVAngle(transI) -SEVAngle(transJ) ).*transG + sin( SEVAngle(transI) -SEVAngle(transJ) ).*transB), ...
varNum,varNum);%线路有功海森矩阵
d2TPij_dTidTj=sparse(transI+busNum,transJ+busNum, ...
Whz(length(mVolt)+2*length(lineI)+1:length(mVolt)+2*length(lineI)+length(transI))...
.*(-SEVolt(transI).*SEVolt(transJ)./transK.*( cos( SEVAngle(transI) -SEVAngle(transJ) ).*transG + sin( SEVAngle(transI) -SEVAngle(transJ) ).*transB) ), ...
varNum,varNum);%线路有功海森矩阵
d2TPij_dTj2=sparse(transJ+busNum,transJ+busNum, ...
Whz(length(mVolt)+2*length(lineI)+1:length(mVolt)+2*length(lineI)+length(transI))...
.*SEVolt(transI).*SEVolt(transJ)./transK.*( cos( SEVAngle(transI) -SEVAngle(transJ) ).*transG + sin( SEVAngle(transI) -SEVAngle(transJ) ).*transB), ...
varNum,varNum);%线路有功海森矩阵
% d2TPij_dTi2=sparse(transI+busNum,transI+busNum, ...
% 1 ...
% .*SEVolt(transI).*SEVolt(transJ)./transK.*( cos( SEVAngle(transI) -SEVAngle(transJ) ).*transG + sin( SEVAngle(transI) -SEVAngle(transJ) ).*transB), ...
% varNum,varNum);%线路有功海森矩阵
% d2TPij_dTidTj=sparse(transI+busNum,transJ+busNum, ...
% 1 ...
% .*(-SEVolt(transI).*SEVolt(transJ)./transK.*( cos( SEVAngle(transI) -SEVAngle(transJ) ).*transG + sin( SEVAngle(transI) -SEVAngle(transJ) ).*transB) ), ...
% varNum,varNum);%线路有功海森矩阵
% d2TPij_dTj2=sparse(transJ+busNum,transJ+busNum, ...
% 1 ...
% .*SEVolt(transI).*SEVolt(transJ)./transK.*( cos( SEVAngle(transI) -SEVAngle(transJ) ).*transG + sin( SEVAngle(transI) -SEVAngle(transJ) ).*transB), ...
% varNum,varNum);%线路有功海森矩阵
% t2=d2TPij_dTi2+d2TPij_dTidTj+d2TPij_dTj2+d2TPij_dTidTj';
% 对电压相角求导
d2TPij_dVidTi=sparse(transI,transI+busNum, ...
Whz(length(mVolt)+length(lineI)*2+1:length(mVolt)+length(lineI)*2+length(transI)) ...
.* (SEVolt(transJ)./transK.* (sin(SEVAngle(transI) -SEVAngle(transJ)).*transG - cos(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
varNum,varNum);%线路有功海森矩阵
d2TPij_dVidTj=sparse(transI,transJ+busNum, ...
Whz(length(mVolt)+length(lineI)*2+1:length(mVolt)+length(lineI)*2+length(transI)) ...
.* (-SEVolt(transJ)./transK.* (sin(SEVAngle(transI) -SEVAngle(transJ)).*transG - cos(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
varNum,varNum);%线路有功海森矩阵
d2TPij_dVjdTi=sparse(transJ,transI+busNum, ...
Whz(length(mVolt)+length(lineI)*2+1:length(mVolt)+length(lineI)*2+length(transI)) ...
.* (SEVolt(transI)./transK.* (sin(SEVAngle(transI) -SEVAngle(transJ)).*transG - cos(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
varNum,varNum);%线路有功海森矩阵
d2TPij_dVjdTj=sparse(transJ,transJ+busNum, ...
Whz(length(mVolt)+length(lineI)*2+1:length(mVolt)+length(lineI)*2+length(transI)) ...
.* (-SEVolt(transI)./transK.* (sin(SEVAngle(transI) -SEVAngle(transJ)).*transG - cos(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
varNum,varNum);%线路有功海森矩阵
% 对电压相角求导
% d2TPij_dVidTi=sparse(transI,transI+busNum, ...
% 1 ...
% .* (SEVolt(transJ)./transK.* (sin(SEVAngle(transI) -SEVAngle(transJ)).*transG - cos(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
% varNum,varNum);%线路有功海森矩阵
% d2TPij_dVidTj=sparse(transI,transJ+busNum, ...
% 1 ...
% .* (-SEVolt(transJ)./transK.* (sin(SEVAngle(transI) -SEVAngle(transJ)).*transG - cos(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
% varNum,varNum);%线路有功海森矩阵
% d2TPij_dVjdTi=sparse(transJ,transI+busNum, ...
% 1 ...
% .* (SEVolt(transI)./transK.* (sin(SEVAngle(transI) -SEVAngle(transJ)).*transG - cos(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
% varNum,varNum);%线路有功海森矩阵
% d2TPij_dVjdTj=sparse(transJ,transJ+busNum, ...
% 1 ...
% .* (-SEVolt(transI)./transK.* (sin(SEVAngle(transI) -SEVAngle(transJ)).*transG - cos(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
% varNum,varNum);%线路有功海森矩阵
% t3=d2TPij_dVidTi+d2TPij_dVidTj+d2TPij_dVjdTi+d2TPij_dVjdTj;
% t=t1+t2+t3+t3';
% 对电压求导
d2TQij_dVi2=sparse(transI,transI, ...
-2*transB./transK./transK.*Whz(length(mVolt)+length(lineI)*2+length(transI)+1:length(mVolt)+length(lineI)*2+length(transI)*2),...
varNum,varNum);%变压器无功海森矩阵
d2TQij_dVidVj=sparse(transI,transJ, ...
Whz(length(mVolt)+length(lineI)*2+length(transI)+1:length(mVolt)+length(lineI)*2+length(transI)*2) ...
.*( -(sin(SEVAngle(transI) -SEVAngle(transJ)).*transG - cos(SEVAngle(transI) -SEVAngle(transJ)).*transB)./transK), ...
varNum,varNum);%变压器无功海森矩阵
d2TQij_dVj2=0;%变压器无功海森矩阵
% 对相角求导
d2TQij_dTi2=sparse(transI+busNum,transI+busNum, ...
Whz(length(mVolt)+length(lineI)*2+length(transI)+1:length(mVolt)+length(lineI)*2+length(transI)*2)...
.*SEVolt(transI).*SEVolt(transJ)./transK.*( sin( SEVAngle(transI) -SEVAngle(transJ) ).*transG - cos( SEVAngle(transI) -SEVAngle(transJ) ).*transB), ...
varNum,varNum);%变压器无功海森矩阵
d2TQij_dTidTj=sparse(transI+busNum,transJ+busNum, ...
Whz(length(mVolt)+length(lineI)*2+length(transI)+1:length(mVolt)+length(lineI)*2+length(transI)*2)...
.*(-SEVolt(transI).*SEVolt(transJ)./transK.*( sin( SEVAngle(transI) -SEVAngle(transJ) ).*transG - cos( SEVAngle(transI) -SEVAngle(transJ) ).*transB)), ...
varNum,varNum);%变压器无功海森矩阵
d2TQij_dTj2=sparse(transJ+busNum,transJ+busNum, ...
Whz(length(mVolt)+length(lineI)*2+length(transI)+1:length(mVolt)+length(lineI)*2+length(transI)*2)...
.*SEVolt(transI).*SEVolt(transJ)./transK.*( sin( SEVAngle(transI) -SEVAngle(transJ) ).*transG - cos( SEVAngle(transI) -SEVAngle(transJ) ).*transB), ...
varNum,varNum);%变压器无功海森矩阵
% 对电压相角求导
d2TQij_dVidTi=sparse(transI,transI+busNum, ...
Whz(length(mVolt)+length(lineI)*2+length(transI)+1:length(mVolt)+length(lineI)*2+length(transI)*2)...
.* (-SEVolt(transJ)./transK.* (cos(SEVAngle(transI) -SEVAngle(transJ)).*transG + sin(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
varNum,varNum);%线路有功海森矩阵
d2TQij_dVidTj=sparse(transI,transJ+busNum, ...
Whz(length(mVolt)+length(lineI)*2+length(transI)+1:length(mVolt)+length(lineI)*2+length(transI)*2)...
.* (SEVolt(transJ)./transK.* (cos(SEVAngle(transI) -SEVAngle(transJ)).*transG + sin(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
varNum,varNum);%线路有功海森矩阵
d2TQij_dVjdTi=sparse(transJ,transI+busNum, ...
Whz(length(mVolt)+length(lineI)*2+length(transI)+1:length(mVolt)+length(lineI)*2+length(transI)*2)...
.* (-SEVolt(transI)./transK.* (cos(SEVAngle(transI) -SEVAngle(transJ)).*transG + sin(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
varNum,varNum);%线路有功海森矩阵
d2TQij_dVjdTj=sparse(transJ,transJ+busNum, ...
Whz(length(mVolt)+length(lineI)*2+length(transI)+1:length(mVolt)+length(lineI)*2+length(transI)*2)...
.* (SEVolt(transI)./transK.* (cos(SEVAngle(transI) -SEVAngle(transJ)).*transG + sin(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
varNum,varNum);%线路有功海森矩阵
%% 形成大的海森矩阵
%海森矩阵对角，即
% [ d2h_dV2 0
% 0 d2h_dT2]
hessian14=d2LPij_dVi2+d2LPij_dVidVj+d2LPij_dVidVj'+d2LPij_dVj2 ...
+d2LPij_dTi2+d2LPij_dTidTj+d2LPij_dTidTj'+d2LPij_dTj2 ...
+d2LQij_dVi2+d2LQij_dVidVj+d2LQij_dVidVj'+d2LQij_dVj2 ...
+d2LQij_dTi2+d2LQij_dTidTj+d2LQij_dTidTj'+d2LQij_dTj2 ...;
+d2TPij_dVi2+d2TPij_dVidVj+d2TPij_dVidVj'+d2TPij_dVj2 ...
+d2TPij_dTi2+d2TPij_dTidTj+d2TPij_dTidTj'+d2TPij_dTj2 ...
+d2TQij_dVi2+d2TQij_dVidVj+d2TQij_dVidVj'+d2TQij_dVj2 ...
+d2TQij_dTi2+d2TQij_dTidTj+d2TQij_dTidTj'+d2TQij_dTj2;
hessian2=d2LPij_dVidTi+d2LPij_dVidTj+d2LPij_dVjdTi+d2LPij_dVjdTj ...
+d2LQij_dVidTi+d2LQij_dVidTj+d2LQij_dVjdTi+d2LQij_dVjdTj ...
+d2TPij_dVidTi+d2TPij_dVidTj+d2TPij_dVjdTi+d2TPij_dVjdTj ...
+d2TQij_dVidTi+d2TQij_dVidTj+d2TQij_dVjdTi+d2TQij_dVjdTj;
hessian3=hessian2';
G=H'*W*H+hessian14+hessian2+hessian3;
%G=H'*W*H;
g=H'*Whz;
% gg=H(5:7,:)'*Whz(5:7);
% ggg=H(busNum++length(lineI)+1:busNum+length(lineI)*2,:)'*H(busNum++length(lineI)+1:busNum+length(lineI)*2,:)+d2LPij_dVi2+d2LPij_dVidVj+d2LPij_dVidVj'+d2LPij_dVj2 ...
% +d2LPij_dTi2+d2LPij_dTidTj+d2LPij_dTidTj'+d2LPij_dTj2 ...
% +d2LPij_dVidTi+d2LPij_dVidTj+d2LPij_dVjdTi+d2LPij_dVjdTj ...
% +[d2LPij_dVidTi+d2LPij_dVidTj+d2LPij_dVjdTi+d2LPij_dVjdTj]';
% ggg=H(busNum++length(lineI)+1:busNum+length(lineI)*2,:)'*H(busNum++length(lineI)+1:busNum+length(lineI)*2,:)+d2LQij_dVi2+d2LQij_dVidVj+d2LQij_dVidVj'+d2LQij_dVj2 ...
% +d2LQij_dTi2+d2LQij_dTidTj+d2LQij_dTidTj'+d2LQij_dTj2 ...
% +d2LQij_dVidTi+d2LQij_dVidTj+d2LQij_dVjdTi+d2LQij_dVjdTj ...
% +[d2LQij_dVidTi+d2LQij_dVidTj+d2LQij_dVjdTi+d2LQij_dVjdTj]';
%
% ggg=H(busNum+length(lineI)*2+length(transI)+1:busNum+length(lineI)*2+length(transI)*2,:)'*H(busNum+length(lineI)*2+length(transI)+1:busNum+length(lineI)*2+length(transI)*2,:)+d2TQij_dVi2+d2TQij_dVidVj+d2TQij_dVidVj'+d2TQij_dVj2 ...
% +d2TQij_dTi2+d2TQij_dTidTj+d2TQij_dTidTj'+d2TQij_dTj2 ...
% +d2TQij_dVidTi+d2TQij_dVidTj+d2TQij_dVjdTi+d2TQij_dVjdTj ...
% +[d2TQij_dVidTi+d2TQij_dVidTj+d2TQij_dVjdTi+d2TQij_dVjdTj]';
% hVolt=H(1:busNum,:)'*H(1:busNum,:);
a=G;
b=-g;
% 平衡节点相角恒定
a(length(mVolt)+Balance,:)=0;
a(:,length(mVolt)+Balance)=0;
a=a+sparse(length(mVolt)+Balance,length(mVolt)+Balance,1,size(a,1),size(a,1));
b(length(mVolt)+Balance)=0;
% 平衡节点电压恒定;
a(Balance,:)=0;
a(:,Balance)=0;
a=a+sparse(Balance,Balance,1,size(a,1),size(a,1));
b(Balance)=0;
dX=a\b;
% a(Balance,:)=0;
% a(:,Balance)=0;
% a=a+sparse(Balance,Balance,1,size(a,1),size(a,1));
% b(Balance)=0;
dX=inv(a)*b;
dXStep=1;
% dXStep=Armijo(z,newwordParameter,W,SEVolt,SEVAngle,dX,g );
maxD=max(abs(dX(1:length(mVolt))))
fprintf('max abs g:%f\n',full(max(abs(g))));
SEVolt=SEVolt+dX(1:length(mVolt))*dXStep;
SEVAngle=SEVAngle+dX(length(mVolt)+1:length(mVolt)*2)*dXStep;
lamda=-dX(length(mVolt)*2+1:end);
% % 计算目标函数下降量
% preObjFun=(z-h)'*W*(z-h);
% oldSEVolt=SEVolt;
% oldSEVAngle=SEVAngle;
% 求更新控制指标 p
% p=(objfun( SEVolt,SEVAngle,W,z,newwordParameter ) - objfun( SEVolt+dX(1:length(mVolt))*dXStep,SEVAngle+dX(length(mVolt)+1:length(mVolt)*2)*dXStep,W,z,newwordParameter ) )/( .5*dX'*(mu*dX-g) );
% if p>0 %接受更新
% % 更新变量
% SEVolt=SEVolt+dX(1:length(mVolt))*dXStep;
% SEVAngle=SEVAngle+dX(length(mVolt)+1:length(mVolt)*2)*dXStep;
% lamda=-dX(length(mVolt)*2+1:end);
% mu=mu*max([1/3,1-(2*p-1)^3]);
% v=2;
% else
% mu=mu*v;
% v=2*v;
% end
SEVAngle=SEVAngle+dX(length(mVolt)+1:end)*dXStep;
Iteration=Iteration+1;
% % 更新后目标函数
% h=[SEVolt;SEBranchP;SEBranchQ;SETransP;SETransQ];
% ojbFunDecrease=preObjFun-(z-h)'*W*(z-h)
% if ojbFunDecrease<1e-3
% mu=100000;
% else
% mu=0;
% end
optimalCondition=-g;
optimalCondition(Balance)=0;
optimalCondition(Balance+length(mVolt))=0;

71
公式/公式.tex
View File

@ -221,6 +221,8 @@ ps.已检验过线路的公式。
\par
线路支路功率二阶导数
有功部分
\newline
海森矩阵对电压
\begin{equation}
\begin{aligned}
\frac{\partial^2 P_{12}}{\partial V_1^2}&=
@ -243,7 +245,7 @@ ps.已检验过线路的公式。
\frac{\partial^2 P_{12}}{\partial V_2^2}&=0
\end{aligned}
\end{equation}
海森矩阵对相角
\begin{equation}
\begin{aligned}
\frac{\partial P_{12}}{\partial \theta_1^2}&=
@ -268,7 +270,41 @@ ps.已检验过线路的公式。
\frac{V_1}{k} V_2[cos(\theta_1 - \theta_2)G_{ij}+sin (\theta_1 - \theta_2)B_{ij}]
\end{aligned}
\end{equation}
海森矩阵对电压相角
\begin{equation}
\begin{aligned}
\frac{\partial P_{12}}{\partial V_1 \partial \theta_1}&=
-\frac{V_2}{k}[-sin(\theta_1 - \theta_2)G_{ij}+cos (\theta_1 - \theta_2)B_{ij}] \\
&=\frac{V_2}{k}[sin(\theta_1 - \theta_2)G_{ij}-cos (\theta_1 - \theta_2)B_{ij}]
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
\frac{\partial P_{12}}{\partial V_1 \partial \theta_2}&=
-\frac{V_2}{k}[sin(\theta_1 - \theta_2)G_{ij}-cos (\theta_1 - \theta_2)B_{ij}]
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
\frac{\partial P_{12}}{\partial V_2 \partial \theta_1}&=
-\frac{V_1}{k}[-sin(\theta_1 - \theta_2)G_{ij}+cos (\theta_1 - \theta_2)B_{ij}] \\
&=\frac{V_1}{k}[sin(\theta_1 - \theta_2)G_{ij}-cos (\theta_1 - \theta_2)B_{ij}]
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
\frac{\partial P_{12}}{\partial V_2 \partial \theta_2}&=
-\frac{V_1}{k}[sin(\theta_1 - \theta_2)G_{ij}-cos (\theta_1 - \theta_2)B_{ij}]
\end{aligned}
\end{equation}
无功部分
\newline
海森矩阵对电压
\begin{equation}
\begin{aligned}
\frac{\partial^2 Q_{12}}{\partial V_1 \partial V_1}&=
@ -288,7 +324,7 @@ ps.已检验过线路的公式。
\frac{\partial^2 Q_{12}}{\partial V_2^2}&=0
\end{aligned}
\end{equation}
海森矩阵对电压
\begin{equation}
\begin{aligned}
\frac{\partial^2 Q_{12}}{\partial \theta_1^2}&=
@ -310,6 +346,37 @@ ps.已检验过线路的公式。
&=\frac{V_1}{k}V_2[sin(\theta_1 - \theta_2)G_{ij}-cos (\theta_1 - \theta_2)B_{ij}]
\end{aligned}
\end{equation}
海森矩阵对电压相角
\begin{equation}
\begin{aligned}
\frac{\partial Q_{12}}{\partial V_1 \partial \theta_1}&=
-\frac{V_2}{k}[cos(\theta_1 - \theta_2)G_{ij}+sin (\theta_1 - \theta_2)B_{ij}]
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
\frac{\partial Q_{12}}{\partial V_1 \partial \theta_2}&=
-\frac{V_2}{k}[-cos(\theta_1 - \theta_2)G_{ij}-sin (\theta_1 - \theta_2)B_{ij}] \\
&=\frac{V_2}{k}[cos(\theta_1 - \theta_2)G_{ij}+sin (\theta_1 - \theta_2)B_{ij}]
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
\frac{\partial Q_{12}}{\partial V_2 \partial \theta_1}&=
-\frac{V_1}{k}[cos(\theta_1 - \theta_2)G_{ij}+sin (\theta_1 - \theta_2)B_{ij}]
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
\frac{\partial Q_{12}}{\partial V_2 \partial \theta_2}&=
-\frac{V_1}{k}[-cos(\theta_1 - \theta_2)G_{ij}-sin (\theta_1 - \theta_2)B_{ij}] \\
&=\frac{V_1}{k}[cos(\theta_1 - \theta_2)G_{ij}+sin (\theta_1 - \theta_2)B_{ij}]
\end{aligned}
\end{equation}
以上公式对线路和没有计及接地支路的变压器适用，只是线路中变比$k$为1.
\end{document}

12
公式/内点法公式.tex
View File

@ -15,7 +15,7 @@ f(x)=[z-h(x)]^T W [z-h(x)]
\end{equation}
最优条件等价为
\begin{equation}
\bigtriangledown f(x)=J^T W [z-h(x)]=0
\bigtriangledown f(x)=J^T W [h(x)-z]=0
\end{equation}
其中$J$是$h(x)$的$Jacobi$矩阵
\newline
@ -55,6 +55,16 @@ f(x)=[z-h(x)]^T W [z-h(x)]
J ^T W J + \tilde{Q}
\end{equation}
其中
$J$具有以下形式
\begin{equation}
J=\left[ %左括号
\begin{array}{cc} %该矩阵一共3列，每一列都居中放置
\frac{\partial f_1}{ \partial x_1} & \frac{\partial f_1}{ \partial x_2}\\ %第一行元素
\frac{\partial f_2}{ \partial x_1} & \frac{\partial f_2}{ \partial x_2}\\ %第二行元素
\end{array}
\right] %右括号
\end{equation}
\begin{equation}
J ^T W J=
(\omega_1 \frac{ \partial f_1 }{\partial x_1} \frac{ \partial f_1 }{\partial x_1}