这里使用的是scipy中的optimize优化函数
具体代码如下：
from scipy import optimize as opt
import numpy as np
def func(X,G):
sum=0.0
for i in np.arange(G.shape[0]):
sum=sum+X[3*i]
return -1*sum
def sigma(X,G,i):
sum=0.0
for j in np.arange(G.shape[0]):
sum=sum+G[i,j]*X[3*j+2]
return sum*X[3*i+2]
def Cons1(G):
cons=()
for i in np.arange(G.shape[0]):
cons=np.append(cons,{'type':'eq','fun':lambda X:X[3*i]-X[3*i+1]-sigma(X,G,i)})
return cons
def Cons2(G,I):
cons=()
for i in np.arange(G.shape[0]):
for j in np.arange(i,G.shape[0]):
cons=np.append(cons,{'type':'ineq','fun':lambda X:I[i,j]-(X[3*i+2]-X[3*j+2])*G[i,j]})
cons=np.append(cons,{'type':'ineq','fun':lambda X:(X[3*i+2]-X[3*j+2])*G[i,j]+I[i,j]})
return cons
def bou(G,P,vn,Pg):
bounds=np.mat(np.zeros([3*G.shape[0],2]))
for i in np.arange(G.shape[0]):
bounds[3*i,:]=Pg[i,:]
bounds[3*i+1,0]=0
bounds[3*i+1,1]=P[0,i]
bounds[3*i+2,0]=0.95*vn
bounds[3*i+2,1]=1.05*vn
return bounds
def x_start(G,P,vn,Pg):
x0=np.array(np.zeros([1,3*G.shape[0]]))
for i in np.arange(G.shape[0]):
x0[0,i]=float((Pg[i,0]+Pg[i,1]))/2
x0[0,i+1]=P[0,i]
x0[0,i+2]=vn
return x0
def opf(G,I,P,vn,Pg):
b=bou(G,P,vn,Pg)
cons=Cons1(G)
cons=np.append(cons,Cons2(G,I))
x=x_start(G,P,vn,Pg)
res=opt.minimize(func, args=(G),x0=x,bounds=b,constraints=cons,method='SLSQP',options={'disp': True})
return res.x
#以下为输入数据：
G=np.mat([[1/0.64,0,-1/0.64,0,0],
[0,1/0.8+1/0.8+1/0.5,-1/0.8,-1/0.8,-1/0.5],
[-1/0.64,-1/0.8,1/0.64+1/0.8+1/1.6,-1/1.6,0],
[0,-1/0.8,-1/1.6,1/1.6+1/0.8+1,-1],
[0,-1/0.5,0,-1,1/0.5+1]])
I=np.mat([[0,0,1000,0,0],
[0,0,1000,1000,1000],
[1000,1000,0,1000,0],
[0,1000,1000,0,1000],
[0,1000,0,1000,0]])
P=np.mat([0,0,1000,800,0.1])
vn=250
Pg=np.mat([[300,400],[1100,1200],[0,0],[0,0],[200,300]])
print(opf(G,I,P,vn,Pg))
