使用WLS进行电力系统状态估计

使用WLS进行电力系统状态估计,加权最小二乘法的电力系统状态估计(http://www.apollocode.net/a/1068.html)

% Power System State Estimation using Weighted Least Square Method..

num = 30; % IEEE - 14 or IEEE - 30 bus system..(for IEEE-14 bus system replace 30 by 14)...

ybus = ybusppg(num); % Get YBus..

zdata = zdatas(num); % Get Measurement data..

bpq = bbusppg(num); % Get B data..

nbus = max(max(zdata(:,4)),max(zdata(:,5))); % Get number of buses..

type = zdata(:,2); % Type of measurement, Vi - 1, Pi - 2, Qi - 3, Pij - 4, Qij - 5, Iij - 6..

z = zdata(:,3); % Measuement values..

fbus = zdata(:,4); % From bus..

tbus = zdata(:,5); % To bus..

Ri = diag(zdata(:,6)); % Measurement Error..

V = ones(nbus,1); % Initialize the bus voltages..

del = zeros(nbus,1); % Initialize the bus angles..

E = [del(2:end); V];   % State Vector..

G = real(ybus);

B = imag(ybus);

vi = find(type == 1); % Index of voltage magnitude measurements..

ppi = find(type == 2); % Index of real power injection measurements..

qi = find(type == 3); % Index of reactive power injection measurements..

pf = find(type == 4); % Index of real powerflow measurements..

qf = find(type == 5); % Index of reactive powerflow measurements..

nvi = length(vi); % Number of Voltage measurements..

npi = length(ppi); % Number of Real Power Injection measurements..

nqi = length(qi); % Number of Reactive Power Injection measurements..

npf = length(pf); % Number of Real Power Flow measurements..

nqf = length(qf); % Number of Reactive Power Flow measurements..

iter = 1;

tol = 5;

while(tol > 1e-4)

    %Measurement Function, h

    h1 = V(fbus(vi),1);

    h2 = zeros(npi,1);

    h3 = zeros(nqi,1);

    h4 = zeros(npf,1);

    h5 = zeros(nqf,1);

    

    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

    

    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

    

    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

    

    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];

    

    % Residue..

    r = z - h;

    

    % Jacobian..

    % 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

    

    % Measurement Jacobian, H..

    H = [H11 H12; H21 H22; H31 H32; H41 H42; H51 H52];

    

    % Gain Matrix, Gm..

    Gm = H'*inv(Ri)*H;

    

    %Objective Function..

    J = sum(inv(Ri)*r.^2);  

    

    % State Vector..

    dE = inv(Gm)*(H'*inv(Ri)*r);

    E = E + dE;

    del(2:end) = E(1:nbus-1);

    V = E(nbus:end);

    iter = iter + 1;

    tol = max(abs(dE));

end

CvE = diag(inv(H'*inv(Ri)*H)); % Covariance matrix..

Del = 180/pi*del;

E2 = [V Del]; % Bus Voltages and angles..

disp('-------- State Estimation ------------------');

disp('--------------------------');

disp('| Bus |    V   |  Angle  | ');

disp('| No  |   pu   |  Degree | ');

disp('--------------------------');

for m = 1:n

    fprintf('%4g', m); fprintf('  %8.4f', V(m)); fprintf('   %8.4f', Del(m)); fprintf('\n');

end

disp('---------------------------------------------');

详情http://www.apollocode.net/a/1068.html

热门文章

暂无图片
编程学习 ·

exe4j详细使用教程(附下载安装链接)

一、exe4j介绍 ​ exe4j是一个帮助你集成Java应用程序到Windows操作环境的java可执行文件生成工具,无论这些应用是用于服务器,还是图形用户界面(GUI)或命令行的应用程序。如果你想在任务管理器中及Windows XP分组的用户友好任务栏…
暂无图片
编程学习 ·

AUTOSAR从入门到精通100讲(126)-浅谈车载充电系统通信方案

01 引言 本文深入研究车载充电系统策略,设计出一套基于电动汽车电池管理系统与车载充电机的CAN通信协议,可供电动汽车设计人员参考借鉴。 02 电动汽车充电系统通讯网络 电动汽车整车控制系统中采用的是CAN总线通信方式,由一个整车内部高速CAN网络、内部低速CAN网络和一个充电…
暂无图片
编程学习 ·

CMake(九):生成器表达式

当运行CMake时,开发人员倾向于认为它是一个简单的步骤,需要读取项目的CMakeLists.txt文件,并生成相关的特定于生成器的项目文件集(例如Visual Studio解决方案和项目文件,Xcode项目,Unix Makefiles或Ninja输入文件)。然…
暂无图片
编程学习 ·

47.第十章 网络协议和管理配置 -- 网络配置(八)

4.3.3 route 命令 路由表管理命令 路由表主要构成: Destination: 目标网络ID,表示可以到达的目标网络ID,0.0.0.0/0 表示所有未知网络,又称为默认路由,优先级最低Genmask:目标网络对应的netmaskIface: 到达对应网络,应该从当前主机哪个网卡发送出来Gateway: 到达非直连的网络,…
暂无图片
编程学习 ·

元宇宙技术基础

请看图: 1、通过AR、VR等交互技术提升游戏的沉浸感 回顾游戏的发展历程,沉浸感的提升一直是技术突破的主要方向。从《愤怒的小鸟》到CSGO,游戏建模方式从2D到3D的提升使游戏中的物体呈现立体感。玩家在游戏中可以只有切换视角,进而提升沉浸…
暂无图片
编程学习 ·

flink的伪分布式搭建

一 flink的伪分布式搭建 1.1 执行架构图 1.Flink程序需要提交给 Job Client2.Job Client将作业提交给 Job Manager3.Job Manager负责协调资源分配和作业执行。 资源分配完成后,任务将提交给相应的 Task Manage。4.Task Manager启动一个线程以开始执行。Task Manage…
暂无图片
编程学习 ·

十进制正整数与二进制字符串的转换(C++)

Function one: //十进制数字转成二进制字符串 string Binary(int x) {string s "";while(x){if(x % 2 0) s 0 s;else s 1 s;x / 2;}return s; } Function two: //二进制字符串变为十进制数字 int Decimal(string s) {int num 0, …
暂无图片
编程学习 ·

[含lw+源码等]微信小程序校园辩论管理平台+后台管理系统[包运行成功]Java毕业设计计算机毕设

项目功能简介: 《微信小程序校园辩论管理平台后台管理系统》该项目含有源码、论文等资料、配套开发软件、软件安装教程、项目发布教程等 本系统包含微信小程序做的辩论管理前台和Java做的后台管理系统: 微信小程序——辩论管理前台涉及技术:WXML 和 WXS…
暂无图片
编程学习 ·

树莓派驱动DHT11温湿度传感器

1,直接使用python库 代码如下 import RPi.GPIO as GPIO import dht11 import time import datetimeGPIO.setwarnings(True) GPIO.setmode(GPIO.BCM)instance dht11.DHT11(pin14)try:while True:result instance.read()if result.is_valid():print(ok)print(&quo…
暂无图片
编程学习 ·

ELK简介

ELK简介 ELK是三个开源软件的缩写,Elasticsearch、Logstash、Kibana。它们都是开源软件。不过现在还新增了一个 Beats,它是一个轻量级的日志收集处理工具(Agent),Beats 占用资源少,适合于在各个服务器上搜集日志后传输给 Logstas…
暂无图片
编程学习 ·

Linux 基础

通常大数据框架都部署在 Linux 服务器上,所以需要具备一定的 Linux 知识。Linux 书籍当中比较著名的是 《鸟哥私房菜》系列,这个系列很全面也很经典。但如果你希望能够快速地入门,这里推荐《Linux 就该这么学》,其网站上有免费的电…
暂无图片
编程学习 ·

Windows2022 无线网卡装不上驱动

想来 Windows2022 和 windows10/11 的驱动应该差不多通用的,但是死活装不上呢? 搜一下,有人提到 “默认安装时‘无线LAN服务’是关闭的,如果需要开启,只需要在“添加角色和功能”中,选择开启“无线LAN服务…
暂无图片
编程学习 ·

【嵌入式面试宝典】版本控制工具Git常用命令总结

目录 创建仓库 查看信息 版本回退 版本检出 远程库 Git 创建仓库 git initgit add <file> 可反复多次使用&#xff0c;添加多个文件git commit -m <message> 查看信息 git status 仓库当前的状态git diff 差异对比git log 历史记录&#xff0c;提交日志--pret…
暂无图片
编程学习 ·

用Postman生成测试报告

newman newman是一款基于nodejs开发的可以运行postman脚本的工具&#xff0c;使用Newman&#xff0c;可以直接从命令运行和测试postman集合。 安装nodejs 下载地址&#xff1a;https://nodejs.org/en/download/ 选择自己系统相对应的版本内容进行下载&#xff0c;然后傻瓜式安…
暂无图片
编程学习 ·

Java面向对象之多态、向上转型和向下转型

文章目录前言一、多态二、引用类型之间的转换Ⅰ.向上转型Ⅱ.向下转型总结前言 今天继续Java面向对象的学习&#xff0c;学习面向对象的第三大特征&#xff1a;多态&#xff0c;了解多态的意义&#xff0c;以及两种引用类型之间的转换&#xff1a;向上转型、向下转型。  希望能…