💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

---

💥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.

🌈4 Matlab代码实现