多端口柔性互联配电网VSC有功无功的电压时序优化【matlab通用版】

程序名称##

多端口柔性互联配电网VSC有功无功电压时序优化【matlab通用版】

程序功能(对象)

1. 网络拓扑适配性：
   1.适用于任意大小的交流电网,支持节点和支路的增删，相关参数修改；
   2.适用于任意大小的直流电网,支持节点和支路的增删，相关参数修改；
2. 分布式电源接入：允许接入风光储等DG的时序出力，且接入数量无上限；
3. VSC的通用性：VSC的端口数需要提前确定，接入交流和直流的位置任意可调；
4. VSC控制策略:控制策略可选Udc-Q/Pac-Uac/Pac-Qac等控制策略，控制策略详见博文：基于交替迭代法的通用型交直流潮流计算方法(matlab版+python版)；在此若考虑不同控制策略优化时间会成倍增加，建议固定损耗率2%.
5. 目标函数：
   F(x) = totloss+detalU (目标函数因人而异)
   totloss:所有交流系统和直流系统、VSC全天损耗综合
   detalU:所有交流节点和直流节点全天电压偏差之和，额定值为1.0p.u.
6. 决策变量：
   1)主控VSC的无功出力(本算例1个)
   2)其余从控VSC的有功和无功出力(本算例2个)
7. 约束条件：潮流约束、电压约束、VSC容量、VSC有功出力约束等
8. 通用型分析：
   1、交流电网、直流电网拓扑允许随意变动；
   2、交流、直流电网均可接入风光储等多个分布式电源；
   3、全天24小时各个光伏和风电出力参数可调，通过VSC互联实现负荷均衡；
   4、IEEE33节点系统33个节点重要程度系数可调，可选择不同的重要程度系数；
9. 计算结果存储：5*24=120个变量存储在individual.mat 方便制图
10. 设置优化前后程序进行对比：优化后电压偏差和网损显著减小；
11. 程序拓展：
    1.可与OTLC/CB、光伏逆变器、储能联合调解无功；
    2.可配合集群划分进行分区电压调控；基于GA+节点法考虑综合性能指标的集群划分【matlab通用版】
    3.可利用多端口实现网络重构进而进行重构优化；

程序函数说明

matlab版：

1. 主程序：main
2. 目标函数子程序：infitness文件
3. 粒子群算法子程序：PSO_2D文件
4. IEEE33节点系统子程序：交流子系统1ac1_case文件、交流子系统2ac2_case文件、交流子系统3ac3_case文件、直流系统dc_case7(本程序未用到，简化直流部分)
5. 牛拉法基波子程序：NLPowerFlow文件
6. 节点导纳子程序：createYbus文件
7. 交流雅克比矩阵子程序：Jacobi文件(n-1+m*n-1+m)
8. 交直流潮流子程序：ACDC_Flow1文件(简化直流潮流)、ACDC_Flow(文件简化计算慢)
   图一为多端口柔性互联配电网VSC有功无功电压时序优化列表
   

程序计算步骤及流程图

以matlab版程序作为讲解。

1 参数初始化

参数初始化，读取网络参数，在此标幺化，无功优化是在基波潮流计算基础之上进行的，参数初始化和基波潮流学习： https://blog.csdn.net/WConstelltion/article/details/123751611.
包括：交流节点负荷，支路阻抗，发电机以及分布式电源参数初始化，直流节点负荷、支路阻抗、VSC各端口控制策略；风光时序出力数据，支路电阻以及分布式电源参数初始。
图二为不同场景下光伏和风电时序出力

图三是不同子系统负荷波动曲线
图四是不同子系统储能出力波动曲线

图五为交直流混联系统，其中采用三端口VSC将三个IEEE33节点系统进行互联。

2主程序计算流程

2.1读取时序参数：(场景需求可以根据自己选择)

%读取典型日实时基准数据中节点负荷参数和分布电源参数
[System1] = xlsread('系统光伏风电以及负荷变化曲线',1);
[System2] = xlsread('系统光伏风电以及负荷变化曲线',2);    %各系统参数表
[System3] = xlsread('系统光伏风电以及负荷变化曲线',3);
total_System = cat(3,System1,System2,System3);  %维度需要保持一致

2.2确定决策变量变化区间：

    PVSC2_range = repmat( [-0.5,0.5] , T , 1 );%正值代表交流侧流入VSC直流
    QVSC2_range = repmat( [-1, 1] , T , 1 );%剩余无功出力系数，不超过最大容量
    PVSC3_range = repmat( [-0.5,0.5] , T , 1 );
    QVSC3_range = repmat( [-1, 1] , T , 1 );
    QVSC1_range = repmat( [-1, 1] , T , 1 );

2.3确定学习因子和惯性权重：

    c1 = c1s+(1-acos(1-i/max_gen)/pi)*(c1e-c1s);
    c2 = c2s+(1-acos(1-i/max_gen)/pi)*(c2e-c2s);
    w = w_max-(w_max-w_min)*i/max_gen;

3 目标函数优化算法计算

%目的：通过算法求解电压偏差和网损最小时的调控策略
3.1依据抽样粒子120维拆分各时段VSC的出力情况，VSC1的有功是要通过潮流计算确定的。

    PVSC2 = individual(1:1*T); Q2max = sqrt(Smax^2 - abs(PVSC2).^2);%无功最大出力范围0.8MVA
    QVSC2 = individual(1*T+1:2*T).*Q2max;
    PVSC3 = individual(2*T+1:3*T); Q3max = sqrt(Smax^2 - abs(PVSC3).^2);%无功最大出力范围0.8MVA
    QVSC3 = individual(3*T+1:4*T).*Q3max; 
    Q1 = individual(4*T+1:5*T);%无功系数

3.2修正各时段负荷、光伏、风电、储能的出力

        %修正交流子系统DG实时出力参数                  
        ac1_data.dg(find(ac1_data.dg(:, 3)==1), 4) = System1(ti , 3)/10^3;%风电修正
        ac1_data.dg(find(ac1_data.dg(:, 3)==2), 4) = System1(ti , 4)/10^3;%光伏修正
        ac1_data.dg(find(ac1_data.dg(:, 3)==3), 4) = -System1(ti , 5)/10^3;%储能作为负荷为正值故需要添加负号      
        ac1_data.dg(find(ac1_data.dg(:, 3)==1), 5) = tan(acos(0.9))*ac1_data.dg(find(ac1_data.dg(:, 3)==1), 4);
        ac1_data.dg(find(ac1_data.dg(:, 3)==2), 5) = tan(acos(0.9))*ac1_data.dg(find(ac1_data.dg(:, 3)==2), 4);
        
        ac2_data.dg(find(ac2_data.dg(:, 3)==1), 4) = System2(ti , 3)/10^3;%风电修正
        ac2_data.dg(find(ac2_data.dg(:, 3)==2), 4) = System2(ti , 4)/10^3;%光伏修正
        ac2_data.dg(find(ac2_data.dg(:, 3)==3), 4) = -System2(ti , 5)/10^3;%储能作为负荷为正值故需要添加负号      
        ac2_data.dg(find(ac2_data.dg(:, 3)==1), 5) = tan(acos(0.9))*ac2_data.dg(find(ac2_data.dg(:, 3)==1), 4);
        ac2_data.dg(find(ac2_data.dg(:, 3)==2), 5) = tan(acos(0.9))*ac2_data.dg(find(ac2_data.dg(:, 3)==2), 4);        
        
        ac3_data.dg(find(ac3_data.dg(:, 3)==1), 4) = System3(ti , 3)/10^3;%风电修正
        ac3_data.dg(find(ac3_data.dg(:, 3)==2), 4) = System3(ti , 4)/10^3;%光伏修正
        ac3_data.dg(find(ac3_data.dg(:, 3)==3), 4) = -System3(ti , 5)/10^3;%储能作为负荷为正值故需要添加负号      
        ac3_data.dg(find(ac3_data.dg(:, 3)==1), 5) = tan(acos(0.9))*ac3_data.dg(find(ac3_data.dg(:, 3)==1), 4);
        ac3_data.dg(find(ac3_data.dg(:, 3)==2), 5) = tan(acos(0.9))*ac3_data.dg(find(ac3_data.dg(:, 3)==2), 4);         

3.3交直流潮流计算：

  [loss1,loss2,loss3,Ss,Uac]=ACDC_Flow(dc_data,ac1_data,ac2_data,ac3_data);

3.4未简化的交直流潮流计算：（详见：[基于交替迭代法的通用型交直流潮流计算方法(matlab版+python版)]），需要反复迭代交直流部分进而确定Pvsc1的大小以及节点电压和网损。简化后的直流部分计算如下所示：(直接通过端口网损率确定)

    Ps = VSC(:,6);Qs = VSC(:,7);Us = VSC(:,8);         %所有换流站注入交流系统的功率  VSC给的是无功系数 
    Pdcvsc = Ps - abs(Ps)*0.02;
    Ps(1) = -sum(Pdcvsc)-abs(sum(Pdcvsc))*0.02;
    Qs(1) = sqrt(Smax^2-(abs(Ps(1)))^2)*Qs(1);
    Ss = Ps(1)+1i*Qs(1);

3.5交流潮流计算各子系统输出参数：

    %牛顿拉夫逊潮流计算，返回节点注入Pi,Qi,节点电压ac_bus,交流迭代次数kac    
    %这里需要进行3个子系统解列求解，概率潮流计算否则不收敛！！！
    [Pac1,Qac1,Bus_V1,kac1]=NLPowerFlow(ac1_baseMVA, ac1_bus,ac1_branch,ac1_gen,ac1_dg);%概率潮流
    Uac1_m = Bus_V1(:,2); Uac1_a = Bus_V1(:,3); 
    [Pac2,Qac2,Bus_V2,kac2]=NLPowerFlow(ac2_baseMVA, ac2_bus,ac2_branch,ac2_gen,ac2_dg);
    Uac2_m = Bus_V2(:,2); Uac2_a = Bus_V2(:,3); 
    [Pac3,Qac3,Bus_V3,kac3]=NLPowerFlow(ac3_baseMVA, ac3_bus,ac3_branch,ac3_gen,ac3_dg);
    Uac3_m = Bus_V3(:,2); Uac3_a = Bus_V3(:,3); 

3.6 计算交直流系统综合损耗,交流子系统损耗，直流系统损耗，换流站损耗：

S1_branch = acbranch_power(ac1_bus,ac1_branch,ac1_baseMVA,Bus_V1,Y1bus);
loss1 = sum(S1_branch(:, 5));
S2_branch = acbranch_power(ac2_bus,ac2_branch,ac2_baseMVA,Bus_V2,Y2bus);
loss2 = sum(S2_branch(:, 5));
S3_branch = acbranch_power(ac3_bus,ac3_branch,ac3_baseMVA,Bus_V3,Y3bus);
loss3 = sum(S3_branch(:, 5));
S_branch = [S1_branch(17,3) S2_branch(17,3) S3_branch(17,3)];%观看传输功率变化

3.7 计算目标函数：(看个人需求，可只要电压偏差)

   totloss = sum(Loss,2);detalU = sum(sum(sum(abs(US-1))));

3.8 主程序绘制VSC的有功无功出力，以及电压网损变化。

程序计算结果分析

1. 1）基于PSO算法迭代曲线如下图所示，未进行优化前目标函数值为99，优化之后目标函数减低值89左右。
   

2）优化前各系统各节点电压、网损变化

3）优化后各子系统节点电压和网损：（网损和电压偏差都是有所减低的，这里需要进一步进行细致对比）

4. matlab版制作不易非无偿，python版有需求可以联系开发，上述仅为博主设置算例，各位硕导博导可以联系博主定制对应的需求。
5. 参考文献：计及源荷不确定性的柔性互联配电系统两阶段电压优化_蔡晖

程序开发

完整版matlab版/python版
lemonyoungman