传统监测点可观测域法的Matlab代码实现

传统监测点可观测域法的Matlab代码实现

clear; close all; clc;

%% 壹、读取算例
example = case39; % 导入matpower6.0中的例子，这里是ieee30节点矩阵
% 使用matpower6.0工具箱读取例子
mpc_example = loadcase(example);

%% 贰、进行稳态潮流计算
% 使用Newtown-Raphson法进行最大误差为1e-10、最大迭代步骤为1000步的交流潮流计算，且不打印输出结果
opt = mpoption('PF_ALG',1,'PF_TOL',1e-10,'PF_MAX_IT',1000,'PF_DC',0,'OUT_ALL',0); % 设置算法类型和迭代条件
fname = 'static power flow for IEEE-39 Power System.txt';
pf = runpf(example,opt,fname); % 进行潮流计算，并将结果打印在txt文件中
% 导出节点导纳矩阵和节点阻抗矩阵
Ybus = cell2mat(struct2cell(load('Ybus.mat'))); % 节点导纳矩阵是稀疏矩阵
Zbus = inv(full(Ybus));% 通过求逆得到节点阻抗矩阵

%% 叁、将计算得到的节点电压数据进行预处理
U_i0 = pf.bus(:,[1,8,9]);% 短路故障前各节点的电压U_i|0|
rad_i0 = (pi/180)*U_i0(:,3);% 短路前故障相角的角度制转换为弧度制

%% 肆、初始化参数
branch = pf.branch;
L = size(branch,1);% 统计线路的数量
Point = 41;% 每条线路平均有41个分段点，也就是将线路平均分成40段
% F = point;
F = L*Point;% 故障点数目等于线路与分段点的乘积
N = size(U_i0,1);
U_th = 0.8;% 电压暂降阈值的标幺值设为0.8
z_f = 0;% 短路点过渡阻抗设置为0

%% 伍、计算故障点f处发生三相短路时各个节点的电压数据
U_f0 = zeros(F,1);% 不同编号的故障点f发生短路故障前的故障点处电压值初始化
CMRA = zeros(F,N);% CMRA矩阵的初始化
Z_if = zeros(F,N);% 不同编号的故障点f发生短路故障时，互阻抗矩阵Z_if的初始化
Z_ff = zeros(F,1);% 不同编号的故障点f发生短路故障时，自阻抗矩阵Z_ff的初始化
U_is = cell(F,1);% 不同编号的故障点f发生短路故障后，各个母线的残压
% 计算故障点和母线之间的互阻抗Z_if以及故障点的自阻抗Z_ff
for line = 1:L
    start = branch(line,1);% 找出线路的起始点
    terminal = branch(line,2);% 找出线路的终止点
    for search = 1:Point
        f = (line-1)*Point+search;% 对离散故障点f进行编号
        lambda = (search-1)/(Point-1);% 计算线路的归一化距离
        U_f0(f) = (1-lambda)*U_i0(start,2)*...
            exp(1j*rad_i0(start))+lambda*U_i0(terminal,2)*...
            exp(1j*rad_i0(terminal));% 计算故障发生前故障点f处的电压值
        Z_if(f,:) = (1-lambda)*Zbus(start,:)+lambda*Zbus(terminal,:);% 计算故障点f与母线i的互阻抗Z_fi
        Z_ff(f) = (1-lambda)^2*Zbus(start,start)+lambda^2*Zbus(terminal,terminal)+...
            2*lambda*(1-lambda)*Zbus(start,terminal)+lambda*(1-lambda)*(branch(line,3)+...
            j*branch(line,4));% 计算故障点f的Z_ff的自阻抗Z_ff
    end
end
% 计算三相短路故障发生后各母线的电压暂降情况
for f = 1:F
    U_is{f} = zeros(N,3);% 对元胞数组的每个元素进行初始化
    U_is{f}(:,1) = U_i0(:,1);% 对母线进行编号
    temp1 = U_i0(:,2).*exp(j*(rad_i0));% temp1为故障前母线电压的复数值
    temp2 = temp1-U_f0(f)*(Z_if(f,:).'/(Z_ff(f)+z_f));% temp2为故障后母线电压的复数值
    U_is{f}(:,2) = abs(temp2);% 故障发生后各母线电压的幅值
    U_is{f}(:,3) = angle(temp2)*180/pi;% 故障发生后各母线电压的相角
    CMRA(f,:) = (U_is{f}(:,2)<=U_th)';
end
clear temp1 temp2;

%% 陆、CMRA优化模型的建立
monitor = @(x) sum(x);% 用于监测点求和，使得监测点数目最少
intcon = 1:N;% 39个节点应全部是整数
reach = ones(F,1);% 应保证所有故障点所有类型的电压暂降源均能被配置的监测点观测到
lb = zeros(N,1);% 下限设为全0
ub = ones(N,1);% 上限设为全1

%% 柒、使用遗传算法(Genetic Algorithm, GA)求解CMRA优化模型
options = gaoptimset('PopulationSize',150);% 设种群规模是150
Iteration_Times = 100;% 使用100次遗传算法进行迭代计算
SolutionSpace = zeros(1,N);% 存放遗传算法的解空间
search_begin = 1;% 初始化粘接处始端变量
search_end = 0;% 初始化粘接处末端变量
% 将遗传算法迭代得到的结果粘接起来
for run = 1:Iteration_Times
    % 使用GA算法算出配置结果
    [x,sum,exitflag,output,population,score] = ga(monitor,N,-CMRA,-reach,[],[],lb,ub,[],intcon,options);
    % 进一步筛选得到满足目标函数和约束条件的所有配置方案
    population = population((score == sum),:);
    % 进行一次粘接
    search_end = search_end+size(population,1);% 移动粘接末端变量
    SolutionSpace(search_begin:search_end,:) = population;% 粘接
    % 一次粘接完成后，准备进行下一次粘接
    search_begin = search_end+1;% 移动下一次粘接需要的始端变量
end
SolutionSpace = delete_rp(SolutionSpace);% 删除解空间中的重复元素
Redundancy = CMRA*SolutionSpace';% 求出各个配置方案的冗余度，验证是否满足约束条件

%% 捌、整合得到监测点优化配置方案
ROW = size(SolutionSpace,1);% 用于对求出的解空间进行按行搜寻
Configuration_Plan = zeros(1,sum);% 预置内存空间
% 整合求出解空间中配置方法对应的监测点配置方案
for row = 1:ROW
    Configuration_Plan(row,:) = find(SolutionSpace(row,:)~=0);
end