文章复现《热电联产系统在区域综合能源系统中的定容选址研究》（matlab程序）

function f=res_pf(x) %x=[5 130] alfa=0.79;%电热比 guanlian=[2 5;8 12;13 15;14 30;11 21];%关联矩阵 xh=guanlian(x(1),1);%热位置 xp=guanlian(x(1),2);%电位置 Pe = x(2); He=Pe/alfa; %% 热网 n=17; % 17节点 m=16; % 16支路 Cp=4.182*10^(-3); %水的比热容，单位：MJ/(kg·℃) % 热力系统负荷节点数据,分别为节点号，节点负荷热功率MW，压力kPa，供热温度，输出温度和回热温度。节点1为平衡节点 ND=[ 1 0 500 100 50 50 2 0 493 100 50 0 3 0.12 414 100 50 50 4 0.12 406 100 50 50 5 0 445 100 50 0 6 0.12 337 100 50 50 7 0.12 345 100 50 50 8 0 395 100 50 0 9 0.12 265 100 50 50 10 0.12 280 100 50 50 11 0 370 100 50 0 12 0.085 337 100 50 50 13 0.085 335 100 50 50 14 0 365 100 50 0 15 0.085 337 100 50 50 16 0.085 335 100 50 50 17 0.085 360 100 50 50]; ND(xh,2)=ND(xh,2)-He*1e-3; % 热力系统管道参数，分别为管道号、始点、终点、长度、直径、传热系数、粗糙度 BD=[1 1 2 70.6 125 0.321 0.4 2 2 3 95.2 32 0.189 0.4 3 2 4 105.1 32 0.189 0.4 4 2 5 61.8 80 0.278 0.4 5 5 6 136 32 0.189 0.4 6 5 7 123.3 32 0.189 0.4 7 5 8 52.1 65 0.236 0.4 8 8 9 161.1 32 0.189 0.4 9 8 10 134.2 32 0.189 0.4 10 8 11 41.7 65 0.236 0.4 11 11 12 136.4 32 0.189 0.4 12 11 13 134.1 32 0.189 0.4 13 11 14 44.9 80 0.278 0.4 14 14 15 136.4 32 0.189 0.4 15 14 16 136.2 32 0.189 0.4 16 14 17 134.2 32 0.189 0.4]; l = BD(:,4); % 长度 m d = BD(:,5)/1000; % 直径 m lmd = BD(:,6); % 传热系数 W/mK rgh = BD(:,7)/1000; % 钢管当量绝对粗糙度m Ta = 10; % 环境温度 10℃ Lh = ND(2:n,2); % 各负荷节点热功率 Ts = ND(:,4); % 供热温度 初值100℃ T0 = ND(2:n,5); % 各负荷节点输出温度 Tr = ND(:,6); % 回热温度 Prs = ND(:,3); % 各节点压力初值 Ao = zeros(n,m); % 构建节点-支路增广关联矩阵A0 ( n x m ) for i=1:n for j=1:m if BD(j,2)==i % 节点i为管道j起点 Ao(i,j)=1; elseif BD(j,3)==i % 节点i为管道j终点 Ao(i,j)=-1; end end end for i=1:m if ND(BD(i,2),3)>ND(BD(i,3),3) D_Vtr(i)=1; % 方向向量D_Vtr else D_Vtr(i)=-1; end end As = Ao(2:n,:); %去掉第一行得到基本关联矩阵 % 构造回路关联矩阵，不存在回路--------------------------------------------- s=zeros(m,m); %阻力系数矩阵 for i=1:m s(i,i)=6.88*10^(-9)*l(i)*((rgh(i))^0.25)/0.9584/((d(i))^5.25); %单位kPa/((m^3/h)^2)，忽略局部阻力当量ld（附件：三通、弯头、阀门等）热水密度0.9584 end deltaP=Ao'*Prs; %管道压差kPa for i=1:m G_branch(i,1)=D_Vtr(i)*sqrt(abs(deltaP(i))/s(i,i)); % 支路流量 单位t/h end for i=1:n-1 G_node(i,1)=3.6*Lh(i)/(Cp*(Ts(i)-Tr(i))); % 节点流量 end F = As*G_branch+G_node; % A*Q+Qq=0 pr = 1e-6; % 误差精度 count = 0; while max(abs(F))>pr %----------------------------水力模型求解---------------------------------- optr = s*diag(abs(G_branch)); % 算子 Jac = As*inv(optr)*As'; Inv_Jac = -inv(Jac); ND(2:n,3) = Inv_Jac*G_node + ND(1,3); for i=1:m if ND(BD(i,2),3)>ND(BD(i,3),3) D_Vtr(i)=1; else D_Vtr(i)=-1; end deltaP=Ao'*ND(:,3); G_branch(i,1)=D_Vtr(i)*sqrt(abs(deltaP(i))/s(i,i)); %--------------------------供热温度------------------------------------ if D_Vtr(i)==1 ND(BD(i,3),4)=(ND(BD(i,2),4)-Ta)*exp(-lmd(i)*l(i)/(3.6*10^6*Cp*abs(G_branch(i))))+Ta; %供热温度 else ND(BD(i,2),4)=(ND(BD(i,3),4)-Ta)/exp(-lmd(i)*l(i)/(3.6*10^6*Cp*abs(G_branch(i))))+Ta; %ND(BD(i,2),4)=(ND(BD(i,3),4)+(exp(-lmd(i)*l(i)/(3.6*10^6*Cp*abs(G_branch(i))))-1)*Ta)/(exp(-lmd(i)*l(i)/(3.6*10^6*Cp*abs(G_branch(i))))); end end %--------------------------回热温度------------------------------------ ND(:,6)=Tr; for i=1:m if D_Vtr(i)==1 ND(BD(i,3),6)=(ND(BD(i,2),6)-Ta)/exp(-lmd(i)*l(i)/(3.6*10^6*Cp*abs(G_branch(i))))+Ta; %ND(BD(i,3),6)=(ND(BD(i,2),6)+(exp(-lmd(i)*l(i)/(3.6*10^6*Cp*abs(G_branch(i))))-1)*Ta)/(exp(-lmd(i)*l(i)/(3.6*10^6*Cp*abs(G_branch(i))))); %回热温度 else ND(BD(i,2),6)=(ND(BD(i,3),6)-Ta)*exp(-lmd(i)*l(i)/(3.6*10^6*Cp*abs(G_branch(i))))+Ta;%回热温度 end end for i=1:n-1 G_node(i,1)=3.6*Lh(i)/(Cp*(ND(i+1,4)-ND(i+1,6))); %单位t/h end F=As*G_branch+G_node; count=count+1; end % 结果 Prs=ND(:,3); %各节点水压kPa G_branch=D_Vtr'.*G_branch; %各管道水流量t/h Ts=ND(:,4); %各节点供热温度 T0=ND(:,5); %各节点输出温度，给定为50℃ Tr=ND(:,6); %各节点回热温度 Ds=Cp*G_branch(1)/3.6*(Ts(1)-Tr(1)); %热源供热功率 Ts; Tr; Ds; %% 电网 % 配电网前推回代潮流计算程序 countnum=0; % 计数 BranchData = [ % 支路矩阵，第一列首节点，第二列末节点，三、四列电阻和电抗 1 2 0.0922 0.0470; 2 3 0.4930 0.2511; 3 4 0.3660 0.1864; 4 5 0.3811 0.1941; 5 6 0.8190 0.7070; 6 7 0.1872 0.6188; 7 8 0.7114 0.2351; 8 9 1.0300 0.7400; 9 10 1.0440 0.7400; 10 11 0.1966 0.0650; 11 12 0.3744 0.1238; 12 13 1.4680 1.1550; 13 14 0.5416 0.7129; 14 15 0.5910 0.5260; 15 16 0.7463 0.5450; 16 17 1.2890 1.7210; 17 18 0.7320 0.5740; 2 19 0.1640 0.1565; 19 20 1.5042 1.3554; 20 21 0.4095 0.4784; 21 22 0.7089 0.9373; 3 23 0.4512 0.3083; 23 24 0.8980 0.7091; 24 25 0.8960 0.7011; 6 26 0.2030 0.1034; 26 27 0.2842 0.1447; 27 28 1.0590 0.9337; 28 29 0.8042 0.7006; 29 30 0.5075 0.2585; 30 31 0.9744 0.9630; 31 32 0.3105 0.3619; 32 33 0.3410 0.5302; ]; % 支路，阻抗 NodeData = [ % 节点矩阵 二、三列注入有功、无功 2 100.00 60.00; 3 90.00 40.00; 4 120.00 80.00; 5 60.00 30.00; 6 60.00 20.00; 7 200.00 100.00; 8 200.00 100.00; 9 60.00 20.00; 10 60.00 20.00; 11 45.00 30.00; 12 60.00 35.00; 13 60.00 35.00; 14 120.00 80.00; 15 60.00 10.00; 16 60.00 20.00; 17 60.00 20.00; 18 90.00 40.00; 19 90.00 40.00; 20 90.00 40.00; 21 90.00 40.00; 22 90.00 40.00; 23 90.00 50.00; 24 420.00 200.00; 25 60 0; 26 60.00 25.00; 27 60.00 25.00; 28 60.00 20.00; 29 120.00 70.00; 30 200.00 600.00; 31 150.00 70.00; 32 210.00 100.00; 33 60.00 40.00; ]; % 节点，负荷 NodeData(xp,2)=NodeData(xp,2)-Pe; UB = 12.66; % 电压基准 kV SB = 10; % 功率基准 MVA ZB = UB^2/SB; % 阻抗基准 ohm BranchData(:,[3,4]) = BranchData(:,[3,4]) / ZB; % 阻抗标幺化 NodeData(:,[2,3]) = NodeData(:,[2,3]) / SB / 1000;% 功率标幺化 NN = 33; % 节点数 A0 = zeros(NN); % A0用于回代过程 for n = 1:NN-1 A0(BranchData(n,1),BranchData(n,2)) = 1; end AssociatedMatrix =0; % 形成关联矩阵 for n=2:NN-1 AssociatedMatrix(n,n)=1;% 把关联矩阵的对角线值变成1 temp=BranchData(n-1,1); AssociatedMatrix(n,1:n-1)=AssociatedMatrix(temp,1:n-1); end A0T = A0'; % 形成A0矩阵的转置，用于前推过程 S = [0;-NodeData(:,2)-sqrt(-1)*NodeData(:,3)]; % 形成 S 节电注入的视在功率 ZL = [0;BranchData(:,3)+sqrt(-1)*BranchData(:,4)]; % 形成 ZL 复阻抗 V = ones(NN,1); % 节点电压向量 V(1) = 1; % 各个节点电压赋初值 IL(NN,1) = -conj(S(NN) / V(NN)); % 最末支路电流赋初值 Delta = 1; % 收敛判据赋初值 TempV = V; % 记忆电压赋初值，用于记忆上次迭代结果 while Delta > 1e-8 % 表示误差精度为1e-8 countnum=countnum+1; IN = conj(S ./ V); % 节点注入电流 for n = 1:NN-1 % 电流回代过程 IL(NN-n) = A0(NN-n,NN-n+1:end)*IL(NN-n+1:end) - IN(NN-n); end for n = 2:NN % 电压前推过程 V(n) = A0T(n,1:n-1)*V(1:n-1) - ZL(n)*IL(n); end Delta = max(abs(V-TempV)); % 更新收敛�