基于动态规划的微电网黑启动优化算例的python实现（原创干货）

  为了发小论文，这俩天写python写了一个微电网黑启动优化算例，由于小论文从投稿到录用还需要几个月，出于保护知识产区的考虑，在这里只分享一个简化后的模型。考虑到很多不是电气专业的朋友，这里给出表示拓扑的简单示意图。

  对于拓扑图中每一个开关 Si，其属性用一个三维向量表示[si,Pi,Ci]。si表示开关对应的开关状态，Pi表示该开关对应节点上负载功率，Ci表示人工打开该开关所需成本。

  本算例优化的目标函数为

这里~表示对位求反，简单说就是0 1互换。

考虑电网的频率稳定，约束函数为闭合每一个开关时，其增加的负载功率导致电网一次调频频率的变化绝对值小于0.2Hz,即：

该算例可以实现在约束条件内求得约束函数最小值对应的开关状态，并求得总成本（如下所示）

接下来是python代码的实现：

busname = ['1_0_0','2_0_0','3_0_0','1_1_0','1_1_1','1_1_2','1_1_3','1_2_0','1_2_1','1_2_2',
           '1_2_3','2_1_0','2_2_0','2_3_0']
numbus = len(busname)
P00 = 9
block = 3

for busi in range(numbus):   #交互式导入Pi
    print('输入Pbus'+str(busi)+'=')
    Pbusi = float(input())
    exec('P'+str(busi)+'=Pbusi')

#交互式导入节点修复花费成本
busi = 0
for busi in range(numbus):
    print('输入Cbus'+str(busi)+'=')
    Cbusi = float(input())
    exec('C'+str(busi)+'=Cbusi')

#Main
num_1bus = 3 #与发电机直接相连的开关数量
S_allpos = 2**(len(busname)-num_1bus) #所有初始状态的可能性总数
for i in range(S_allpos):
    S_bin = to14(i)
    S_bin = list(S_bin)
    s1_0_0 = float(S_bin[0])
    s2_0_0 = float(S_bin[1])
    s3_0_0 = float(S_bin[2])
    s1_1_0 = float(S_bin[3])
    s1_1_1 = float(S_bin[4])
    s1_1_2 = float(S_bin[5])
    s1_1_3 = float(S_bin[6])
    s1_2_0 = float(S_bin[7])
    s1_2_1 = float(S_bin[8])
    s1_2_2 = float(S_bin[9])
    s1_2_3 = float(S_bin[10])
    s2_1_0 = float(S_bin[11])
    s2_2_0 = float(S_bin[12])
    s2_3_0 = float(S_bin[13])
    if aaa_S_O(i) == 0:
        exec('M'+str(i)+'=9999')
    else:
        Mi = (C0*int(not(s1_0_0))+C1*int(not(s2_0_0))+
             C2*int(not(s3_0_0))+C3*int(not(s1_1_0))+
             C4*int(not(s1_1_1))+C5*int(not(s1_1_2))+
             C6*int(not(s1_1_3))+C7*int(not(s1_2_0))+
             C8*int(not(s1_2_1))+C9*int(not(s1_2_2))+
             C10*int(not(s1_2_3))+C11*int(not(s2_1_0))+
             C12*int(not(s2_2_0))+C13*int(not(s2_3_0)))
        exec('M'+str(i)+'=Mi')
temp2 = 0
M = 9999
for i in range(S_allpos):
    if eval('M'+str(i)) < eval('M'+str(temp2)):
        temp2 = i
print(to14(temp2))
print(eval('M'+str(temp2)))

#计算有几个开关需要被合上
def count_0(strx):
    count = 0
    for i in strx:
        if i == '0':
            count = count+1
    returnary = []
    i = 0
    count1 = 0
    
    for j in strx:
        count1 = count1 +1
        if j == '0':
            returnary.append(count1-1)                         
    return [count,returnary]

def to14(input1):
#返回第input1种故障状态转化为的14位二进制数
    temp1 = bin(input1)[2:]
    if len(temp1) < 14 :
        for i in range(14-len(temp1)):
            temp1 = '0'+temp1
    return temp1

def Pf(input1):
    returnnum = input1*0.33
    return returnnum 

def aaa_S_O(S):
    S_bin = to14(S)
    S_bin = list(S_bin)
    [count_bus,ary1] = count_0(S_bin) #有几个需要合上的开关及位置
    temp1 = 0 
    s1_0_0 = float(S_bin[0])
    s2_0_0 = float(S_bin[1])
    s3_0_0 = float(S_bin[2])
    s1_1_0 = float(S_bin[3])
    s1_1_1 = float(S_bin[4])
    s1_1_2 = float(S_bin[5])
    s1_1_3 = float(S_bin[6])
    s1_2_0 = float(S_bin[7])
    s1_2_1 = float(S_bin[8])
    s1_2_2 = float(S_bin[9])
    s1_2_3 = float(S_bin[10])
    s2_1_0 = float(S_bin[11])
    s2_2_0 = float(S_bin[12])
    s2_3_0 = float(S_bin[13])
    Pij = (s1_1_1*s1_1_0*s1_0_0*P4+s1_1_2*s1_1_0*s1_0_0*P5
           +s1_1_3*s1_1_0*s1_0_0*P6+s1_2_1*s1_2_0*s1_0_0*P8
           +s1_2_2*s1_2_0*s1_0_0*P9+s1_2_3*s1_2_0*s1_0_0*P10
           +s2_1_0*s2_0_0*P11+s2_2_0*s2_0_0*P12+s2_3_0*s2_0_0*P13
           +s3_0_0*P2)+P00 
    returnnum = 1
    count_bus1 = count_bus
    for j in range(count_bus): 
        temp1 = 0
        for i in range(count_bus1):         
            if caldp(S_bin,temp1)>caldp(S_bin,i):
                temp1 = i
        if caldp(S_bin,temp1)>Pf(Pij):
            returnnum = 0
        
        count_bus1 = count_bus1 -1
        Pij =Pij + caldp(S_bin,temp1)
        S_bin[ary1[temp1]] = 1
        s1_0_0 = float(S_bin[0])
        s2_0_0 = float(S_bin[1])
        s3_0_0 = float(S_bin[2])
        s1_1_0 = float(S_bin[3])
        s1_1_1 = float(S_bin[4])
        s1_1_2 = float(S_bin[5])
        s1_1_3 = float(S_bin[6])
        s1_2_0 = float(S_bin[7])
        s1_2_1 = float(S_bin[8])
        s1_2_2 = float(S_bin[9])
        s1_2_3 = float(S_bin[10])
        s2_1_0 = float(S_bin[11])
        s2_2_0 = float(S_bin[12])
        s2_3_0 = float(S_bin[13])  
        del(ary1[temp1])                       
    return returnnum

def caldp(input1,input2):
#计算在input1状态时input2开关合上后的deltaP
    S_bin = input1
    P00 = 9
    S_bin = list(S_bin)
    s1_0_0 = float(S_bin[0])
    s2_0_0 = float(S_bin[1])
    s3_0_0 = float(S_bin[2])
    s1_1_0 = float(S_bin[3])
    s1_1_1 = float(S_bin[4])
    s1_1_2 = float(S_bin[5])
    s1_1_3 = float(S_bin[6])
    s1_2_0 = float(S_bin[7])
    s1_2_1 = float(S_bin[8])
    s1_2_2 = float(S_bin[9])
    s1_2_3 = float(S_bin[10])
    s2_1_0 = float(S_bin[11])
    s2_2_0 = float(S_bin[12])
    s2_3_0 = float(S_bin[13])
    P_old = (s1_1_1*s1_1_0*s1_0_0*P4+s1_1_2*s1_1_0*s1_0_0*P5
           +s1_1_3*s1_1_0*s1_0_0*P6+s1_2_1*s1_2_0*s1_0_0*P8
           +s1_2_2*s1_2_0*s1_0_0*P9+s1_2_3*s1_2_0*s1_0_0*P10
           +s2_1_0*s2_0_0*P11+s2_2_0*s2_0_0*P12+s2_3_0*s2_0_0*P13
           +s3_0_0*P2)+P00 
    S_bin[input2] = float(1)
    s1_0_0 = float(S_bin[0])
    s2_0_0 = float(S_bin[1])
    s3_0_0 = float(S_bin[2])
    s1_1_0 = float(S_bin[3])
    s1_1_1 = float(S_bin[4])
    s1_1_2 = float(S_bin[5])
    s1_1_3 = float(S_bin[6])
    s1_2_0 = float(S_bin[7])
    s1_2_1 = float(S_bin[8])
    s1_2_2 = float(S_bin[9])
    s1_2_3 = float(S_bin[10])
    s2_1_0 = float(S_bin[11])
    s2_2_0 = float(S_bin[12])
    s2_3_0 = float(S_bin[13])
    P_new = (s1_1_1*s1_1_0*s1_0_0*P4+s1_1_2*s1_1_0*s1_0_0*P5
           +s1_1_3*s1_1_0*s1_0_0*P6+s1_2_1*s1_2_0*s1_0_0*P8
           +s1_2_2*s1_2_0*s1_0_0*P9+s1_2_3*s1_2_0*s1_0_0*P10
           +s2_1_0*s2_0_0*P11+s2_2_0*s2_0_0*P12+s2_3_0*s2_0_0*P13
           +s3_0_0*P2)+P00 
    dp = P_new - P_old
    return dp