💥1 概述
当前世界资源与环境问题日益突出,对能源系统进行合理的安全评估是保障其经济、稳定运行的重要基础。评估结果不仅可以实时反馈系统的运行情况,还能为其优化调度提供指导。然而,随着多种能源系统形成的综合能源体系规模不断增大、耦合日益紧密、拓扑结构日趋复杂,传统的评估方法多用于能流形式单一、负荷模式固定的场景,已不足以对复杂的综合能源系统供能安全性给出准确评估。
📚2 运行结果
部分代码：
clear; clc; tic; %%% 设置基准值 Hload_B=1; %非线性，不要动了 Hm_B=1; % T_B=1; CP_B=1; HK_B=1; G_B=1; Gpi_B=1; Kg_B=1; %% 迭代参数设置 Elimit=1*10^-5; Glimit=10; Hlimit=1*10^-1; Tmax=100; %% CHP机组参数 % cm=823200; cm=755528; ne=0.0015; %% 电网参数输入 mpc_E = loadcase('case_E.m'); E_nb = length(mpc_E.bus); E_n_PQ = sum(mpc_E.bus(:,2) == 1); % PQ 节点个数 P=zeros(1,E_nb); %P,Q为原始数据，Pi,Qi为计算结果 Q=zeros(1,E_nb); U=ones(1,E_nb); %电压初始值由此确定 theta=zeros(1,E_nb); P(mpc_E.bus(:,1)) = - mpc_E.bus(:,3); P(mpc_E.gen(:,1)) = mpc_E.gen(:,2); Q(mpc_E.bus(:,1)) = - mpc_E.bus(:,4); U(mpc_E.gen(:,1)) = mpc_E.gen(:,6); theta(mpc_E.bus(mpc_E.bus(:,2) == 3,1)) = mpc_E.bus(mpc_E.bus(:,2) == 3,9); % theta = mpc_E.bus(,9); P = P/mpc_E.baseMVA; Q = Q/mpc_E.baseMVA; Y = makeYbus(mpc_E.baseMVA,mpc_E.bus,mpc_E.branch); G=real(Y); %电导矩阵 B=imag(Y); %电纳矩阵 %% 热网参数输入 mpc_test=case_H(); H_nb = length(mpc_test.bus); H_nl = length(mpc_test.branch); hkk=0.2;%热传导系数 Cp=4200;%热媒的比热容 Cp1=4200; % 原始算例 % hload=[50000,50000,50000,50000,100000,50000,10000,50000,80000,100000,50000,50000,1000000]; % 增加热功率 % index = (Gpipe(:,5)~=3); % temp = sum(index); % Ag = sparse(Gpipe(index,1),1:temp,-1,G_nb,temp) + sparse(Gpipe(index,2),1:temp,1,G_nb,temp); %输入平衡节点气压 单位bar 负荷节点默认为0 Gp1=[5.000 5.000]/Gpi_B; %天然气系统数据处理 zg = sum(Gpipe(:,5) == 3); %计算管道常数 % Kr = 70840/Kg_B./sqrt(Gpipe(:,3))'; % disp('管道摩擦系数：'); % disp(Kr); %计算迭代初始值，这部分再商量，先假设初值为下述： Gpi=[4.4973 4.4839 4.4397 4.4394 4.4323]/Gpi_B; Gpi=[Gp1,Gpi]; Gpi2=Gpi.^2; % disp('Gpi2:'); % disp(Gpi2); %计算气网不平衡量 [Gload,Gpi2,dpi2,Ag,Gpipe,Ag2,dGQ,dG] = GUnbalanced(0,Gpi2,Ag,Gpipe,G_nl,zg,0,Gload,gas); % 计算热网不平衡量 [Cs,bs,Cr,br,hk,hT0source] = H_un(mh,lh,sh,hkk,Hpipe,HTs,Hm,hT0,HTrload,Cp1); [dP_H,dp_H,dTs,dTr] = bupinghengliang(Cp,As1,Hm,hT0,Bh,HK,Cs,HTsload,bs,Cr,HTrload,br,Hload); %整合为一整个不平衡量 dEGH=[dP';dQ';dP_H;dp_H;dTs;dTr;dG]; % disp('整个的不平衡量：'); % disp(dEGH); %% 牛拉法大循环开始 for round=1:Tmax fprintf('第%d次迭代\n',round); %%%%%%%%%%%耦合设备转换 hT0source=40/T_B; % PG=Cp*(Hm(1)+Hm(4))*(HTssource-hT0source); PG=Cp*Hm(1)*(HTssource-hT0source); P(30)=(PG/cm); Gload(3)=PG/(cm*ne); % 计算雅克比矩阵 % 求解Jee J=Jacobi( E_nb,E_n_PQ,U,theta,B,G,Pi,Qi ); % J = -makeJac();PQ 顺序不同%$% % disp('Jee:'); % disp(J); % 求解Jhh [HJhh] = HJacobi(nh,mh,lh,Cp,HTsload,HTs,hT0,As1,Bh,HK,Hm,Cs,Cr,Hpipe,hk); % 求解Jgg for i=1:G_nl - zg D(i,i)=dGQ(i)/(2*dpi2(i)); end % disp(D); Jgg=-Ag2*D*Ag2'; % disp('Jgg：'); % disp(Jgg); % Jgh(3,1)=Cp*Hm(4)*(HTssource-hT0source)/(cm*ne); % Jgh(3,4)=Cp*Hm(1)*(HTssource-hT0source)/(cm*ne); Jgh(3,1)=Cp*(HTssource-hT0source)/(cm*ne); % disp('Jgh'); % disp(Jgh); % 求解Jge Jge=zeros(5,E_n_PQ+E_nb-1); % 求解Jhe Jhe=zeros(mh+2*lh,E_n_PQ+E_nb-1); % 求解Jhg Jhg=zeros(mh+2*lh,5); % 合成整个的雅克比矩阵 JEGH=[J,Jeh,Jeg; Jhe,HJhh,Jhg; Jge,Jgh,Jgg]; % disp('JEGH'); % disp(JEGH); con(round) = cond(JEGH); con1=max(con); % 求解修正量 [dtheta,dU,dHm,dHTsload,dHTrload,dppp] = EGHcorrect(E_nb,E_n_PQ,mh,lh,U,dEGH,JEGH,5); % 修正状态量 % 修正电网状态量 U=U+dU; theta=theta+dtheta; % disp('节点电压幅值：'); % disp(U); % disp('节点电压相角：'); % disp(rad2deg(theta)); % 修正热网状态量 Hm=Hm+dHm; HTsload=HTsload+dHTsload; HTrload=HTrload+dHTrload; % disp('修正后的管道流量：'); % disp(Hm); % disp('修正后的节点供热温度：'); % disp(HTsload); % disp('修正后的节点回热温度：'); % disp(HTrload); % 修正气网状态量 Gpi2=Gpi2+dppp; % disp('修正后的节点气压值'); % disp(Gpi2); % 判断方向 % [As1,Hpipe,Hm] = panduan(mh,Hpipe,As1,Hm); % 求解不平衡量 % 计算电网不平衡量 [dP,dQ,Pi,Qi]=EUnbalanced( E_nb,E_n_PQ,P,Q,U,G,B,theta ); [ref, pv, pq] = bustypes(mpc_E.bus, mpc_E.gen); V = U' .* exp(theta' * j); Vm = U'; mpopt = mpoption; mis = V .* conj(Y * V) - makeSbus(mpc_E.baseMVA, mpc_E.bus, mpc_E.gen, mpopt, Vm); F = [ real(mis(1:38));%pv;pq;30： imag(mis(pq)) ]; %计算气网不平衡量 [Gload,Gpi2,dpi2,Ag,Gpipe,Ag2,dGQ,dG] = GUnbalanced(0,Gpi2,Ag,Gpipe,G_nl,zg,0,Gload,gas); % 计算热网不平衡量 [Cs,bs,Cr,br,hk,hT0source] = H_un(mh,lh,sh,hkk,Hpipe,HTs,Hm,hT0,HTrload,Cp1); [dP_H,dp_H,dTs,dTr] = bupinghengliang(Cp,As1,Hm,hT0,Bh,HK,Cs,HTsload,bs,Cr,HTrload,br,Hload); %整合为一整个不平衡量 dEGH=[dP';dQ';dP_H;dp_H;dTs;dTr;dG]; % disp('整个的不平衡量：'); % disp(dEGH); % 判断收敛 if (max(abs(dP))<Elimit && max(abs(dQ))<Elimit ) && max(abs(dP_H))<Hlimit && max(abs(dp_H))<Hlimit && max(abs(dTs))<Hlimit && max(abs(dTr))<Hlimit && max(abs(dG))<Glimit break; end end %% 迭代结束，判断收敛，输出结果 if (max(abs(dP))<Elimit && max(abs(dQ))<Elimit ) && max(abs(dP_H))<Hlimit && max(abs(dp_H))<Hlimit && max(abs(dTs))<Hlimit && max(abs(dTr))<Hlimit && max(abs(dG))<Glimit disp('计算收敛'); else disp('计算不收敛'); end disp('************最终结果*************'); fprintf('迭代总次数：%d\n', round); disp('***电网计算结果：') disp('节点电压幅值：'); disp(U); disp('节点电压相角：'); disp(rad2deg(theta)); disp('有功计算结果：'); disp(Pi); disp('无功计算结果：'); disp(Qi); disp('***热网计算结果：'); disp('负荷节点供热温度:'); disp(HTsload); disp('负荷节点回热温度:'); disp(HTrload); disp('热源回热温度:'); disp(hT0source); disp('热网管道流量：'); disp(Hm); disp('***气网计算结果：'); disp('天然气网络各节点压力：'); Gpi=sqrt(Gpi2); disp(Gpi); disp('天然气网络各管道流量：'); disp(dGQ); disp('天然气网络负荷：'); disp(Gload); toc;
🎉3 参考文献
部分理论来源于网络，如有侵权请联系删除。
[1]董洪辛. 综合能源系统供能安全评估及其应用[D].大连理工大学,2022.DOI:10.26991/d.cnki.gdllu.2022.001547.
