Distflow_OPF.zip_Distflow_OPF_distflow 程序_wind power flow_热电_热电潮

%======================================================================== % Function: IEEE 33 system Distflow_OPF % Author: Cao Yang % Date: 2018/01/30 % State: Success %======================================================================== clear clc %% 配电网参数 % 调度周期总时段 NT = 24; T = 1; % 平衡节点电压 Vsl = 1; %% 基准值 Sbase = 100; % MVA Vbase= 30;% kV Zbase = Vbase^2/Sbase;% o Ibase = Sbase/(sqrt(3)*Vbase);% kA %% 节点数据 % 节点编号(1)|节点类型(2)|有功负荷(3)|无功负荷(4)|节点电导(5)|节点电纳(6) busdata = [ 1 3 0 0 0 0; % MW MVar 2 1 0.500 0.460 0 0; 3 1 0.490 0.440 0 0; 4 1 0.520 0.580 0 0; 5 1 0.370 0.430 0 0; 6 1 0.360 0.420 0 0; 7 1 0.500 0.500 0 0; 8 1 0.500 0.500 0 0; 9 1 0.560 0.420 0 0; 10 1 0.460 0.420 0 0; 11 1 0.315 0.530 0 0; 12 1 0.460 0.435 0 0; 13 1 0.440 0.435 0 0; 14 1 0.420 0.380 0 0; 15 1 0.460 0.410 0 0; 16 1 0.560 0.520 0 0; 17 1 0.460 0.420 0 0; 18 1 0.410 0.440 0 0; 19 1 0.590 0.540 0 0; 20 1 0.490 0.440 0 0; 21 1 0.490 0.440 0 0; 22 1 0.490 0.440 0 0; 23 1 0.590 0.450 0 0; 24 1 0.420 0.500 0 0; 25 1 0.420 0.400 0 0; 26 1 0.460 0.525 0 0; 27 1 0.560 0.425 0 0; 28 1 0.450 0.420 0 0; 29 1 0.420 0.570 0 0; 30 1 0.320 0.500 0 0; 31 1 0.550 0.470 0 0; 32 1 0.410 0.500 0 0; 33 1 0.460 0.440 0 0; ]; % 电负荷曲线 Electric_curve = [ 0.45 0.44 0.43 0.45 0.46 0.57 0.65 0.72 0.85 0.93 0.95 1.00 0.98 0.94 0.95 0.95 0.95 0.95 0.91 0.83 0.78 0.78 0.56 0.51 ]; N_bus = size(busdata,1); % 线路节点数 % 有功负荷|无功负荷 P_d = bsxfun(@times,busdata(:,3),Electric_curve'); % MW//用0.45乘以所有行元素形成的是在第一个计时点每个Bus的负荷 即电负荷曲线是每个节点在该时刻满发的百分比 Q_d = bsxfun(@times,busdata(:,4),Electric_curve'); % Mvar % 化为标幺值 P_d = P_d/Sbase; Q_d = Q_d/Sbase; U_bus_min = 0.95; U_bus_max = 1.05; g_bus = busdata(:,5) * Zbase; b_bus = busdata(:,6) * Zbase;%电导和电纳 %% 发电机数据 % 机组位置(1)|最大无功(2)|最小无功(3)|最大有功(4)|最小有功(5)|成本系数A(6)|成本系数B(7)|成本系数C(8) gendata = [ 1 10 -10 10 0.1 1.8697 25.6438 31.67; % MW 3 10 -10 10 0.1 1.8297 26.2438 31.67; 7 10 -10 10 0.1 1.8973 26.1438 31.67; ]; N_gen = size(gendata,1);%发电机数 Pg_max = gendata(:,4)/Sbase; % 标幺值 Pg_min = gendata(:,5)/Sbase; % 标幺值 Qg_max = gendata(:,2)/Sbase; % 标幺值 Qg_min = gendata(:,3)/Sbase; % 标幺值 a_gen = gendata(:,6); % 成本系数A b_gen = gendata(:,7); % 成本系数B c_gen = gendata(:,8); % 成本系数C Index_gen = gendata(:,1); % 发电机节点编号 %分段线性化近似//采用分段线性化近似以后目标函数就不会出现非线性的了 N_piece = 5; % 5段近似 Pk_gen = zeros(N_gen,N_piece+1); Ck_gen = zeros(N_gen,N_piece+1); for i = 1:N_gen Pk_gen(i,:) = linspace(Pg_min(i),Pg_max(i),N_piece+1);%他把发电机的输出功率分成了六段，CKgen是开到每一功率时的成本的行向量 Ck_gen(i,:) = a_gen(i) .* (Pk_gen(i,:) * Sbase).^2 + b_gen(i) .* (Pk_gen(i,:) * Sbase) + c_gen(i);%所有档位对应的发电成本 end %% CHP供热机组 % 电网中机组位置(1)|热网中机组位置(2)|气-电效率(3)|气-热效率(4)|电出力最小值(5)|电出力最大值(6)|热出力最小值(7)|热出力最大值(8) CHPdata = [ 11 15 0.32 0.46 0.1 3 0 3 ]; N_CHP = size(CHPdata,1); ef_CHP_E = CHPdata(:,3); % CHP的电效率 ef_CHP_H = CHPdata(:,4); % CHP的热效率 P_CHP_min = CHPdata(:,5)/Sbase; % CHP的电出力最小值（标幺值） P_CHP_max = CHPdata(:,6)/Sbase; % CHP的电出力的最大值（标幺值） Q_CHP_min = CHPdata(:,7); % CHP的热出力的最小值 Q_CHP_max = CHPdata(:,8); % CHP的热出力的最大值 %% 风电机组数据 % 机组位置(1)|出力预测（2-25） Wind_gendata = [ % MW 16 0.9846 0.9846 0.7692 0.6000 1.0000 0.8231 0.8462 0.8077 0.7692 0.3077 0.1538 0.1923 0.2308 0.1923 0.2077 0.5385 0.8077 0.5385 0.4615 0.6923 0.7692 0.8077 0.8462 0.8462; % MW 29 0.9846 0.9846 0.7692 0.6000 1.0000 0.8231 0.8462 0.8077 0.7692 0.3077 0.1538 0.1923 0.2308 0.1923 0.2077 0.5385 0.8077 0.5385 0.4615 0.6923 0.7692 0.8077 0.8462 0.8462; ]; N_windgen = size(Wind_gendata,1); Pwg_max = Wind_gendata(:,2:(NT+1))/Sbase; % 标幺值/出力的最大值矩阵 Pwg_min = zeros(N_windgen,NT)/Sbase; % 标幺值/出力的最小值矩阵 %% 线路数据 % 线路编号(1)|起始节点(2)|终止节点(3)|线路电阻(4)|线路电抗(5)|有功潮流上限(6)|有功潮流下限(7)|无功潮流上限(8)|无功潮流下限(9) branchdata = [ 1 1 2 0.0922 0.0470 9.9 0 9.9 0; 2 2 3 0.2930 0.2512 9.9 0 9.9 0; 3 3 4 0.3661 0.1864 9.9 0 9.9 0; 4 4 5 0.3811 0.1941 9.9 0 9.9 0; 5 5 6 0.8190 0.7070 9.9 0 9.9 0; 6 6 7 0.1872 0.6188 9.9 0 9.9 0; 7 7 8 0.7115 0.2351 9.9 0 9.9 0; 8 8 9 1.0299 0.7400 9.9 0 9.9 0; 9 9 10 1.0440 0.7400 9.9 0 9.9 0; 10 10 11 0.1967 0.0651 9.9 0 9.9 0; 11 11 12 0.3744 0.1298 9.9 0 9.9 0; 12 12 13 1.4680 1.1549 9.9 0 9.9 0; 13 13 14 0.5416 0.7129 9.9 0 9.9 0; 14 14 15 0.5909 0.5260 9.9 0 9.9 0; 15 15 16 0.7462 0.5449 9.9 0 9.9 0; 16 16 17 1.2889 1.7210 9.9 0 9.9 0; 17 17 18 0.7320 0.5739 9.9 0 9.9 0; 18 2 19 0.1640 0.1565 9.9 0 9.9 0; 19 19 20 1.5042 1.3555 9.9 0 9.9 0; 20 20 21 0.4095 0.4784 9.9 0 9.9 0; 21 21 22 0.7089 0.9373 9.9 0 9.9 0; 22 3 23 0.4512 0.3084 9.9 0 9.9 0; 23 23 24 0.8980 0.7091 9.9 0 9.9 0; 24 24 25 0.8959 0.7071 9.9 0 9.9 0; 25 6 26 0.2031 0.1034 9.9 0 9.9 0; 26 26 27 0.2842 0.1447 9.9 0 9.9 0; 27 27 28 1.0589 0.9338 9.9 0 9.9 0; 28 28 29 0.8043 0.7006 9.9 0 9.9 0; 29 29 30 0.5074 0.2585 9.9 0 9.9 0; 30 30 31 0.9745 0.9629 9.9 0 9.9 0; 31 31 32 0.3105 0.3619 9.9 0 9.9 0; 32 32 33 0.3411 0.5302 9.9 0 9.9 0; ]; N_line = size(branchdata,1); % 系统线路数 % 电阻和电抗(p.u.) r_line = branchdata(:,4)/Zbase; x_line = branchdata(:,5)/Zbase; Pmax_line = branchdata(:,6)/Sbase; Pmin_line = branchdata(:,7)/Sbase; Qmax_line = branchdata(:,8)/Sbase; Qmin_line = branchdata(:,9)/Sbase; i_bus_line = branchdata(:,2);%from bus j_bus_line = branchdata(:,3);%to bus % 建立每个节点上下游的索引集合 for i = 1: N_bus temp1 = find(i_bus_line == i); % 线路首节点为i/// 搜索出的是首节点的对应线路，即首节点的下游线路 temp2 = find(j_bus_line == i); % 线路末节点为i////搜索出的是末节点的对应线路，即末节点的上游线路 Up_line(i,1:length(temp2)) = temp2; % 节点i的上游线路集合///因为upline共有33行，所以每一行对应一个节点，故每一行都是对应temp2的末节点，所以temp2是i的上游线路 Down_line(i,1:length(temp1)) = temp1; % 节点i的下游线路集合//下游线路同理 Up_bus(i,1:length(temp2)) = branchdata(temp2,2); % 节点i的上游节点集合//i是temp2的末节点，且线路temp2的首节点为 branchdata(temp2,2)，所以首节点是末节点的上游节点 Down_bus(i,1:length(temp1)) = branchdata(temp1,3); % 节点i的下游节点集合///用length的原因是 i可能是多条线路的首末节点 若i只是一条线路的首末节点，那么其他列系统会自动填成0 end %% 定义变量 P_line = sdpvar(N_line,NT,'full'); % 线路有功 Q_line = sdpvar(N_line,NT,'full'); % 线路无功 U_bus= sdpvar(N_bus,NT,'full'); % 节点电压幅值 P_g = sdpvar(N_gen,NT,'full'); % 各节点发电机注入有功 Q_g = sdpvar(N_gen,NT,'full'); % 各节点发电机注入无功 P_wg = sdpvar(N_windgen,NT,'full'); % 各节点风电注入有功 P_CHP = sdpvar(N_CHP,NT,'full'); % CHP的电出力 lam_gen = sdpvar(N_piece+1,NT,N_gen,'full'); % 发电机成本近似线性化的系�