程序代码WLS.zip_wls_加权最小二乘_状态估计_程序代码WLS

% WLS电力系统状态估计 num = 14; % IEEE - 14 or IEEE - 30 bus system..... ybus = ybusppg(num); % 得到 YBus.. zdata = zdatas(num); % 得到 测量数据.. bpq = bbusppg(num); % 得到 B 数据.. nbus = max(max(zdata(:,4)),max(zdata(:,5))); % 节点数.. type = zdata(:,2); % 测量类型, Vi - 1, Pi - 2, Qi - 3, Pij - 4, Qij - 5, Iij - 6.. z = zdata(:,3); % 测量值.. %v1=1; %v2=100; %a=0.05; %M=41; %N=1; %x=gaussmix(M,N,a,v1,v2); %z=z+x; z = z + 0.036*randn(length(z),1); fbus = zdata(:,4); % 起始节点.. tbus = zdata(:,5); % 终止节点.. Ri = diag(zdata(:,6)); % 测量误差.. V = ones(nbus,1); % 初始化总线电压.. del = zeros(nbus,1); % 初始化角度.. E = [del(2:end); V]; % 状态向量.. G = real(ybus); B = imag(ybus); tru=tr(num); vi = find(type == 1); % 电压幅值测量.. ppi = find(type == 2); % 实际注入功率.. qi = find(type == 3); % 注入无功功率.. pf = find(type == 4); % 实际潮流.. qf = find(type == 5); % 无功潮流.. nvi = length(vi); % 电压测量数.. npi = length(ppi); % 实际注入功率测量数.. nqi = length(qi); % 无功功率注入测量数. npf = length(pf); % 实际潮流测量数.. nqf = length(qf); % 无功潮流测量数.. iter = 1; tol = 5; %true tv=tru(:,2); tdel=tru(:,3); %for mcc ker = 2;%ker=s V_mcc = ones(nbus,1); % 初始化总线电压.. del_mcc = zeros(nbus,1); % 初始化相角.. E_mcc = [del_mcc(2:end); V_mcc]; % 状态向量.. while(tol > 1e-6) %测算值, h h1 = V(fbus(vi),1);%母线i电压幅值 h2 = zeros(npi,1); h3 = zeros(nqi,1); h4 = zeros(npf,1); h5 = zeros(nqf,1); %h2 母线i有功注入功率 for i = 1:npi m = fbus(ppi(i)); for k = 1:nbus h2(i) = h2(i) + V(m)*V(k)*(G(m,k)*cos(del(m)-del(k)) + B(m,k)*sin(del(m)-del(k))); end end %h3 母线i无功注入功率 for i = 1:nqi m = fbus(qi(i)); for k = 1:nbus h3(i) = h3(i) + V(m)*V(k)*(G(m,k)*sin(del(m)-del(k)) - B(m,k)*cos(del(m)-del(k))); end end %h4 支路mn有功潮流 for i = 1:npf m = fbus(pf(i)); n = tbus(pf(i)); h4(i) = -V(m)^2*G(m,n) - V(m)*V(n)*(-G(m,n)*cos(del(m)-del(n)) - B(m,n)*sin(del(m)-del(n))); end %h5 支路mn无功潮流 for i = 1:nqf m = fbus(qf(i)); n = tbus(qf(i)); h5(i) = -V(m)^2*(-B(m,n)+bpq(m,n)) - V(m)*V(n)*(-G(m,n)*sin(del(m)-del(n)) + B(m,n)*cos(del(m)-del(n))); end h = [h1; h2; h3; h4; h5]; % [zz1 zz2] = size(h); %z1 = h+0.05*randn(zz1,zz2); % 残量.. %r = z - z1; r = z - h; % 雅克比矩阵.. % H11 - Derivative of V with respect to angles.. All Zeros H11 = zeros(nvi,nbus-1); % H12 - Derivative of V with respect to V.. H12 = zeros(nvi,nbus); for k = 1:nvi for n = 1:nbus if n == k H12(k,n) = 1; end end end % H21 - Derivative of Real Power Injections with Angles.. H21 = zeros(npi,nbus-1); for i = 1:npi m = fbus(ppi(i)); for k = 1:(nbus-1) if k+1 == m for n = 1:nbus H21(i,k) = H21(i,k) + V(m)* V(n)*(-G(m,n)*sin(del(m)-del(n)) + B(m,n)*cos(del(m)-del(n))); end H21(i,k) = H21(i,k) - V(m)^2*B(m,m); else H21(i,k) = V(m)* V(k+1)*(G(m,k+1)*sin(del(m)-del(k+1)) - B(m,k+1)*cos(del(m)-del(k+1))); end end end % H22 - Derivative of Real Power Injections with V.. H22 = zeros(npi,nbus); for i = 1:npi m = fbus(ppi(i)); for k = 1:(nbus) if k == m for n = 1:nbus H22(i,k) = H22(i,k) + V(n)*(G(m,n)*cos(del(m)-del(n)) + B(m,n)*sin(del(m)-del(n))); end H22(i,k) = H22(i,k) + V(m)*G(m,m); else H22(i,k) = V(m)*(G(m,k)*cos(del(m)-del(k)) + B(m,k)*sin(del(m)-del(k))); end end end % H31 - Derivative of Reactive Power Injections with Angles.. H31 = zeros(nqi,nbus-1); for i = 1:nqi m = fbus(qi(i)); for k = 1:(nbus-1) if k+1 == m for n = 1:nbus H31(i,k) = H31(i,k) + V(m)* V(n)*(G(m,n)*cos(del(m)-del(n)) + B(m,n)*sin(del(m)-del(n))); end H31(i,k) = H31(i,k) - V(m)^2*G(m,m); else H31(i,k) = V(m)* V(k+1)*(-G(m,k+1)*cos(del(m)-del(k+1)) - B(m,k+1)*sin(del(m)-del(k+1))); end end end % H32 - Derivative of Reactive Power Injections with V.. H32 = zeros(nqi,nbus); for i = 1:nqi m = fbus(qi(i)); for k = 1:(nbus) if k == m for n = 1:nbus H32(i,k) = H32(i,k) + V(n)*(G(m,n)*sin(del(m)-del(n)) - B(m,n)*cos(del(m)-del(n))); end H32(i,k) = H32(i,k) - V(m)*B(m,m); else H32(i,k) = V(m)*(G(m,k)*sin(del(m)-del(k)) - B(m,k)*cos(del(m)-del(k))); end end end % H41 - Derivative of Real Power Flows with Angles.. H41 = zeros(npf,nbus-1); for i = 1:npf m = fbus(pf(i)); n = tbus(pf(i)); for k = 1:(nbus-1) if k+1 == m H41(i,k) = V(m)* V(n)*(-G(m,n)*sin(del(m)-del(n)) + B(m,n)*cos(del(m)-del(n))); else if k+1 == n H41(i,k) = -V(m)* V(n)*(-G(m,n)*sin(del(m)-del(n)) + B(m,n)*cos(del(m)-del(n))); else H41(i,k) = 0; end end end end % H42 - Derivative of Real Power Flows with V.. H42 = zeros(npf,nbus); for i = 1:npf m = fbus(pf(i)); n = tbus(pf(i)); for k = 1:nbus if k == m H42(i,k) = -V(n)*(-G(m,n)*cos(del(m)-del(n)) - B(m,n)*sin(del(m)-del(n))) - 2*G(m,n)*V(m); else if k == n H42(i,k) = -V(m)*(-G(m,n)*cos(del(m)-del(n)) - B(m,n)*sin(del(m)-del(n))); else H42(i,k) = 0; end end end end % H51 - Derivative of Reactive Power Flows with Angles.. H51 = zeros(nqf,nbus-1); for i = 1:nqf m = fbus(qf(i)); n = tbus(qf(i)); for k = 1:(nbus-1) if k+1 == m H51(i,k) = -V(m)* V(n)*(-G(m,n)*cos(del(m)-del(n)) - B(m,n)*sin(del(m)-del(n))); else if k+1 == n H51(i,k) = V(m)* V(n)*(-G(m,n)*cos(del(m)-del(n)) - B(m,n)*sin(del(m)-del(n))); else H51(i,k) = 0; end end end end % H52 - Derivative of Reactive Power Flows with V.. H52 = zeros(nqf,nbus); for i = 1:nqf m = fbus(qf(i)); n = tbus(qf(i)); for k = 1:nbus if k == m H52(i,k) = -V(n)*(-G(m,n)*sin(del(m)-del(n)) + B(m,n)*cos(del(m)-del(n))) - 2*V(m)*(-B(m,n)+ bpq(m,n)); else if k == n H52(i,k) = -V(m)*(-G(m,n)*sin(del(m)-del(n)) + B(m,n)*cos(del(m)-del(n))); else H52(i,k) = 0; end end end end % 测量雅可比.. H = [H11 H12; H21 H22; H31 H32; H41 H42; H51 H52]; % 增益矩阵, Gm.. Gm = H'*inv(Ri)*H; %for wls t1 = 1/(ker^2); W1 = t1*exp(-r.^2/(2*ker^2)); W = diag(W1); Wt = diag(r.^2/(ker^2)); Gmcc = H'*inv(Ri)*W*(Wt-eye(size(Wt,1)))*H; dE_mcc = -inv(Gmcc)*(H'*inv(Ri)*W*r); E_mcc = E_mcc + dE_mcc; del_mcc(2:end) = E_mcc(1:nbus-1); V_mcc = E_mcc(nbus:end); J_mcc = sum(inv(Ri)*exp(-r.^2/(2*ker^2))); %目标函数.. J = sum(inv(Ri)*r.^2); %状态向量.. dE = inv(Gm)*(H'*inv(Ri)*r); E = E + dE;