|
|
★
小木虫: 金币+0.5, 给个红包，谢谢回帖
最近正在编写的一段有关一个小型电网的拓扑搜索计算程序，以下我贴出我编写的拓扑搜索分支程序，请各位大牛多多赐教！
subroutine PowerSystemParameterGetValue()
use PowerSystemConstant
use PowerSystemParameter
use Comp_TP
implicit none
logical alive
! PhysicalNodeNum=8
! GenNum=2
! Tran2Num=1
! BreakerNum=2
! IsolatorNum=2
! GroundCapacitanceNum=1
! LoadNum=2
! BusNum=4
! LineNum=3
inquire(file="Data/PowerSystemParameter.txt",exist=alive)
if (alive)then
open(unit=FileId,file='Data/PowerSystemParameter.txt')
read(FileId,*)PhysicalNodeNum,GenNum,Tran2Num,BreakerNum,IsolatorNum,GroundCapacitanceNum,LoadNum,BusNum,LineNum
write(*,*)" Read PowerSystemParameter.txt is OK!"
close(FileId,status='keep')
else
write(*,*)" PowerSystemParameter.txt doesn't exist."
endif
endsubroutine
subroutine Com_TPAllocate()
use PowerSystemParameter
use Comp_TP
implicit none
if (BreakerNum/=0)then
allocate(TP_Breakers(BreakerNum))
endif
if (IsolatorNum/=0)then
allocate(TP_Isolators(IsolatorNum))
endif
if (BusNum/=0)then
allocate(TP_Buss(BusNum))
endif
if (LoadNum/=0)then
allocate(TP_Loads(LoadNum))
endif
if (GroundCapacitanceNum/=0)then
allocate(TP_GroundCapacitances(GroundCapacitanceNum))
endif
if (LineNum/=0)then
allocate(TP_Lines(LineNum))
endif
if (GenNum/=0)then
allocate(TP_Gens(GenNum))
endif
if (Tran2Num/=0)then
allocate(TP_Tran2s(Tran2Num))
endif
endsubroutine
!=========================================
! 拓扑分析元件数据初始化
!=========================================
subroutine NodeInit()
!拓扑搜索赋值
use PowerSystemConstant
use PowerSystemParameter
use Comp_TP
implicit none
integer::i
logical alive
inquire(file="Data/Breaker.txt",exist=alive)
if (alive)then
open(unit=FileId,file='Data/Breaker.txt')
do i=1,BreakerNum
read(FileId,*)TP_Breakers(i).P1,TP_Breakers(i).P2,TP_Breakers(i).isoff
enddo
close(FileId,status='keep')
write(*,*)" Read Breaker.txt is OK!"
!write(*,*)TP_Breakers(BreakerNum)
else
write(*,*)" Breaker.txt doesn't exist."
endif
inquire(file="Data/Isolator.txt",exist=alive)
if (alive)then
open(unit=FileId,file='Data/Isolator.txt')
do i=1,IsolatorNum
read(FileId,*)TP_Isolators(i).P1,TP_Isolators(i).P2,TP_Isolators(i).isoff
enddo
close(FileId,status='keep')
write(*,*)" Read Isolator.txt is OK!"
! write(*,*)TP_Isolators(IsolatorNum)
else
write(*,*)" Isolator.txt doesn't exist."
endif
inquire(file="Data/Bus.txt",exist=alive)
if (alive)then
open(unit=FileId,file='Data/Bus.txt')
do i=1,BusNum
read(FileId,*)TP_Buss(i).Enable,TP_Buss(i).P1
enddo
close(FileId,status='keep')
write(*,*)" Read Bus.txt is OK!"
! write(*,*)TP_Buss(BusNum)
else
write(*,*)" Bus.txt doesn't exist."
endif
inquire(file="Data/Load.txt",exist=alive)
if (alive)then
open(unit=FileId,file='Data/Load.txt')
do i=1,LoadNum
read(FileId,*)TP_Loads(i).Enable,TP_Loads(i).P1
enddo
close(FileId,status='keep')
write(*,*)" Read Load.txt is OK!"
!write(*,*)TP_Loads(LoadNum)
else
write(*,*)" Load.txt doesn't exist."
endif
inquire(file="Data/GroundCapacitance.txt",exist=alive)
if (alive)then
open(unit=FileId,file='Data/GroundCapacitance.txt')
do i=1,GroundCapacitanceNum
read(FileId,*)TP_GroundCapacitances(i).Enable,TP_GroundCapacitances(i).P1
enddo
close(FileId,status='keep')
write(*,*)" Read GroundCapacitance.txt is OK!"
!write(*,*)TP_GroundCapacitances(GroundCapacitanceNum)
else
write(*,*)" GroundCapacitance.txt doesn't exist."
endif
inquire(file="Data/Generator.txt",exist=alive)
if (alive)then
open(unit=FileId,file='Data/Generator.txt')
do i=1,GenNum
read(FileId,*)TP_Gens(i).Enable,TP_Gens(i).P1,TP_Gens(i).NodeType
enddo
close(FileId,status='keep')
write(*,*)" Read Generator.txt is OK!"
!write(*,*)TP_Gens(GenNum)
else
write(*,*)" Generator.txt doesn't exist."
endif
inquire(file="Data/Tran2.txt",exist=alive)
if (alive)then
open(unit=FileId,file='Data/Tran2.txt')
do i=1,Tran2Num
read(FileId,*)TP_Tran2s(i).Enable,TP_Tran2s(i).P1,TP_Tran2s(i).P2
enddo
close(FileId,status='keep')
write(*,*)" Read Tran2.txt is OK!"
!write(*,*)TP_Tran2s(Tran2Num)
else
write(*,*)" Tran2.txt doesn't exist."
endif
inquire(file="Data/Line.txt",exist=alive)
if (alive)then
open(unit=FileId,file='Data/Line.txt')
do i=1,LineNum
read(FileId,*)TP_Lines(i).Enable,TP_Lines(i).P1,TP_Lines(i).P2
enddo
close(FileId,status='keep')
write(*,*)" Read Line.txt is OK!"
! write(*,*)TP_Lines(LineNum)
else
write(*,*)" Line.txt doesn't exist."
endif
endsubroutine
!=====================================================
! 全网逻辑节点形成函数
!=====================================================
subroutine GetLogicalNode()
use PowerSystemParameter
use TopologyAanlysisVars
use Comp_TP
implicit none
integer(kind=4) i,j,max,min
!动态开辟物理节点到逻辑节点数组（位于TopologyAanlysisVars模块）内存空间
allocate(PhyicalToLogicalArray(PhysicalNodeNum))
!数组初始化
do i=1,PhysicalNodeNum
PhyicalToLogicalArray(i)=i
enddo
!隔离开关状态量处理——隔离开关融合
if (IsolatorNum/=0)then
do i=1,IsolatorNum
if (TP_Isolators(i).isoff==1)then
if(PhyicalToLogicalArray(TP_Isolators(i).P1)>PhyicalToLogicalArray(TP_Isolators(i).P2))then
max=PhyicalToLogicalArray(TP_Isolators(i).P1)
min=PhyicalToLogicalArray(TP_Isolators(i).P2)
else
max=PhyicalToLogicalArray(TP_Isolators(i).P2)
min=PhyicalToLogicalArray(TP_Isolators(i).P1)
endif
PhyicalToLogicalArray(TP_Isolators(i).P1)=min
PhyicalToLogicalArray(TP_Isolators(i).P2)=min
if (max/=min)then
do j=1,PhysicalNodeNum
if(PhyicalToLogicalArray(j)==max)then
PhyicalToLogicalArray(j)=min
elseif(PhyicalToLogicalArray(j)>max)then
PhyicalToLogicalArray(j)=PhyicalToLogicalArray(j)-1
endif
enddo
endif
endif
enddo
endif
!断路器状态量处理——断路器融合
if (BreakerNum/=0)then
do i=1,BreakerNum
if (TP_Breakers(i).isoff==1)then
if(PhyicalToLogicalArray(TP_Breakers(i).P1)>PhyicalToLogicalArray(TP_Breakers(i).P2))then
max=PhyicalToLogicalArray(TP_Breakers(i).P1)
min=PhyicalToLogicalArray(TP_Breakers(i).P2)
else
max=PhyicalToLogicalArray(TP_Breakers(i).P2)
min=PhyicalToLogicalArray(TP_Breakers(i).P1)
endif
PhyicalToLogicalArray(TP_Breakers(i).P1)=min
PhyicalToLogicalArray(TP_Breakers(i).P2)=min
if (max/=min)then
do j=1,PhysicalNodeNum
if(PhyicalToLogicalArray(j)==max)then
PhyicalToLogicalArray(j)=min
elseif(PhyicalToLogicalArray(j)>max)then
PhyicalToLogicalArray(j)=PhyicalToLogicalArray(j)-1
endif
enddo
endif
endif
enddo
endif
!查找网络逻辑节点最大值
LogicalNodeNum= PhyicalToLogicalArray(1)
do i=1,PhysicalNodeNum
if( PhyicalToLogicalArray(i)>LogicalNodeNum)then
LogicalNodeNum=PhyicalToLogicalArray(i)
endif
enddo
write(*,*)"逻辑节点数组为:"
write(*,*)PhyicalToLogicalArray
write(*,*)"逻辑节点最大值为:"
write(*,*)LogicalNodeNum
endsubroutine
!================================================================
! 网络电气岛处理函数
!================================================================
subroutine GetIsland()
use PowerSystemParameter
use TopologyAanlysisVars
use Comp_TP
implicit none
integer(kind=4)::i,j,max,min,SaveIsLandNum
integer(kind=4)::iTemp,jTemp
if(LogicalNodeNum/=0)then
allocate(IsLandArray(LogicalNodeNum))
endif
if(LogicalNodeNum/=0)then
do i=1,LogicalNodeNum
IsLandArray(i)=i
enddo
endif
!输电线路电气岛融合
if(LineNum/=0)then
do i=1,LineNum
if(IsLandArray(PhyicalToLogicalArray(TP_Lines(i).P1))>IsLandArray(PhyicalToLogicalArray(TP_Lines(i).P2)))then
max=IsLandArray(PhyicalToLogicalArray(TP_Lines(i).P1))
min=IsLandArray(PhyicalToLogicalArray(TP_Lines(i).P2))
else
max=IsLandArray(PhyicalToLogicalArray(TP_Lines(i).P2))
min=IsLandArray(PhyicalToLogicalArray(TP_Lines(i).P1))
endif
IsLandArray(PhyicalToLogicalArray(TP_Lines(i).P1))=min
IsLandArray(PhyicalToLogicalArray(TP_Lines(i).P2))=min
if(max/=min)then
do j=1,LogicalNodeNum
if(IsLandArray(j)==max)then
IsLandArray(j)=min
elseif(IsLandArray(j)>max)then
IsLandArray(j)=IsLandArray(j)-1
endif
enddo
endif
enddo
endif
!双绕组变压器电气岛融合
if(Tran2Num/=0)then
do i=1,Tran2Num
if(IsLandArray(PhyicalToLogicalArray(TP_Tran2s(i).P1))>IsLandArray(PhyicalToLogicalArray(TP_Tran2s(i).P2)))then
max=IsLandArray(PhyicalToLogicalArray(TP_Tran2s(i).P1))
min=IsLandArray(PhyicalToLogicalArray(TP_Tran2s(i).P2))
else
max=IsLandArray(PhyicalToLogicalArray(TP_Tran2s(i).P2))
min=IsLandArray(PhyicalToLogicalArray(TP_Tran2s(i).P1))
endif
IsLandArray(PhyicalToLogicalArray(TP_Tran2s(i).P1))=min
IsLandArray(PhyicalToLogicalArray(TP_Tran2s(i).P2))=min
if(max/=min)then
do j=1,LogicalNodeNum
if(IsLandArray(j)==max)then
IsLandArray(j)=min
elseif(IsLandArray(j)>max)then
IsLandArray(j)=IsLandArray(j)-1
endif
enddo
endif
enddo
endif
!===========保留平衡节点所在的电气岛，其他岛不要==================
SaveIsLandNum=0
if(GenNum/=0)then
do i=1,GenNum
if((TP_Gens(i).Enable==1).and.(TP_Gens(i).NodeType==2))then
SaveIsLandNum=IsLandArray(PhyicalToLogicalArray(TP_Gens(i).P1))
! write(*,*)'平衡节点所在岛号'
! write(*,*)SaveIsLandNum
exit
endif
enddo
if(SaveIsLandNum==0)then
write(*,*) "未设置平衡节点！请先设置平衡节点再运行程序"
endif
endif
if(LogicalNodeNum/=0)then
do i=1,LogicalNodeNum
if(SaveIsLandNum==IsLandArray(i))then
IsLandArray(i)=1
else
IsLandArray(i)=0
endif
enddo
endif
!================统计岛上逻辑节点数目==============
NodeNum=0
if(LogicalNodeNum/=0)then
do i=1,LogicalNodeNum
if(IsLandArray(i)==1)then
NodeNum=NodeNum+1
endif
enddo
endif
!===================逻辑节点再排序================
if(NodeNum/=0)then
allocate(IsLandToNodeArray(NodeNum))
allocate(NodeTypeArray(NodeNum))
j=0
do i=1,LogicalNodeNum
if(IsLandArray(i)==1)then
j=j+1
IsLandToNodeArray(j)=i
NodeTypeArray(j)=0 !节点类型初始化，都为PQ节点0，pv=1，slack=2
endif
enddo
endif
if(NodeNum/=0)then
do i=1,NodeNum
do j=1,GenNum
if((TP_Gens(j).NodeType==2).and.(PhyicalToLogicalArray(TP_Gens(j).P1)==IsLandToNodeArray(i)))then
NodeTypeArray(i)=2
elseif((TP_Gens(j).NodeType==1).and.(PhyicalToLogicalArray(TP_Gens(j).P1)==IsLandToNodeArray(i)))then
NodeTypeArray(i)=1
endif
enddo
enddo
endif
!===============节点排序===========
if(NodeNum/=0)then
!==========slack节点排序到最后==============
do j=1,NodeNum
if(NodeTypeArray(j)==2)then
jTemp=NodeTypeArray(j)
NodeTypeArray(j)= NodeTypeArray(NodeNum)
NodeTypeArray(NodeNum)=jTemp
jTemp=IsLandToNodeArray(j)
IsLandToNodeArray(j)=IsLandToNodeArray(NodeNum)
IsLandToNodeArray(NodeNum)=jTemp
endif
enddo
!=====================pv节点排序================
do i=1,NodeNum-1
do j=NodeNum-1,1,-1
if(i
if((NodeTypeArray(i)==1).and.(NodeTypeArray(j)==0))then
iTemp=NodeTypeArray(i)
NodeTypeArray(i)=NodeTypeArray(j)
NodeTypeArray(j)=iTemp
iTemp=IsLandToNodeArray(i)
IsLandToNodeArray(i)=IsLandToNodeArray(j)
IsLandToNodeArray(j)=iTemp
endif
endif
enddo
enddo
! i=1
! do j=1,NodeNum-1
! if(NodeTypeArray(j)==1)then
! jTemp=NodeTypeArray(j)
! NodeTypeArray(j)=NodeTypeArray(NodeNum-i)
! NodeTypeArray(NodeNum-i)=jTemp
! jTemp=IsLandToNodeArray(j)
! IsLandToNodeArray(j)=IsLandToNodeArray(NodeNum-i)
! IsLandToNodeArray(NodeNum-i)=jTemp
! i=i+1
! endif
! enddo
! enddo
endif
!============排序后的物理节点和电气岛上的逻辑节点对应关系======
if(PhysicalNodeNum/=0)then
allocate(PhyicalToNodeArray(PhysicalNodeNum))
do i=1,PhysicalNodeNum
PhyicalToNodeArray(i)=0
enddo
do i=1,NodeNum
do j=1,PhysicalNodeNum
if (IsLandToNodeArray(i)==PhyicalToLogicalArray(j))then
PhyicalToNodeArray(j)=i
endif
enddo
enddo
endif
write(*,*)"电气岛处理结果如下"
write(*,*)IsLandArray
write(*,*)"岛上的逻辑节点编号为："
write(*,*)IsLandToNodeArray
write(*,*)"岛上的逻辑节点编号对应的节点类型为："
write(*,*)NodeTypeArray
write(*,*)"岛上的逻辑节点数目为："
write(*,*)NodeNum
write(*,*)"拓扑搜索的最终成果为："
write(*,*)PhyicalToNodeArray
endsubroutine |
|
回复此楼
送鲜花一朵 
搞定了嘛？回复晚了，抱歉 
