我正在尝试使用GEKKO优化电源系统。具体来说,使用MPC来测试IEEE 14 bus测试用例。
[系统包括14条总线,模型由状态变量theta
和omega
(分别为发电机的功率角和旋转频率)和代数方程组成。
代数方程式对应于DC近似值,即潮流方程的线性近似值。在w^ref
以下的等式中,目标频率为60 Hz。 P
是每条总线上的有功功率矢量。 B
是系统的电纳矩阵。 u
是与发电机的机械功率相对应的控制输入。
目标是通过操纵机械功率w^ref
来使每个总线的0
和功率角接近u
。
我得到的错误是(下面的代码):
in dc_opf
m.solve(disp=True,debug=True) File
"/.local/lib/python2.7/site-packages/gekko/gekko.py", line 1957, in solve
self._build_model() File
"/.local/lib/python2.7/site-packages/gekko/gk_write_files.py", line 33, in _build_model
if not (parameter.VALUE==None): File
"/.local/lib/python2.7/site-packages/gekko/gk_operators.py", line 25, in __len__
return len(self.value) File
"/.local/lib/python2.7/site-packages/gekko/gk_operators.py", line 144, in __len__
return len(self.value)
TypeError: object of type 'int' has no len()
我的问题是编码哪里出错?
我有两个功能dc_opf()
和dc_mats(mat,mode)
。前者是进行优化的地方。后者是用于填充P
和B
矩阵的辅助函数。
我的代码是:
from gekko import GEKKO
import numpy as np
def dc_opf():
m = GEKKO(remote=False)
omega_ref = m.Param(value=60.) #m.Array(m.FV,(14,1))
omega_hi = m.Param(value=61.)
omega_lo = m.Param(value=59.)
H = m.Array(m.FV,(14,1))
Hs = [5.15, 6.54, 6.54, 0., 0., 5.06, 0., 5.06,0.,0.,0.,0.,0.,0.] #Moment of inertia
for i in range(14):
H[i,0].value= Hs[i]
P = m.Array(m.FV,(14,1))
P = dc_mats(P, 'Pow_full')
theta = m.Array(m.SV,(14,1))
u = m.Array(m.CV,(14,1))
for i in range(14):
u[i,0].STATUS = 1
omega = m.Array(m.SV,(14,1))
B = m.Array(m.FV,(14,14))
B = dc_mats(B, 'B_full')
# Soft constraints
oH = m.CV(value=0)
oL = m.CV(value=0)
oH.SPHI=0; oH.WSPHI=100; oH.WSPLO=0 ; oH.STATUS = 1
oL.SPLO=0; oL.WSPHI=0 ; oL.WSPLO=100; oL.STATUS = 1
m.Equations([oH==omega-omega_hi,oL==omega-omega_lo])
m.Equations([theta[i,0].dt() == omega-omega_ref for i in range(14)])
m.Equations([omega[i,0].dt() == (u-P)/(2.0*H) for i in range(14)])
m.Equation(P == B*theta)
m.Minimize((theta) + (omega-omega_ref) + (u-P))
m.options.IMODE = 6
m.solve(disp=True,debug=True)
def dc_mats(mat,mode):
ppc = {"version": '2'}
ppc["baseMVA"] = 100.0 # system MVA base
ppc['branch'] = np.array([
[1, 2, 0.01938, 0.05917, 0.0528, 9900, 0, 0, 0, 0, 1, -360, 360],
[1, 5, 0.05403, 0.22304, 0.0492, 9900, 0, 0, 0, 0, 1, -360, 360],
[2, 3, 0.04699, 0.19797, 0.0438, 9900, 0, 0, 0, 0, 1, -360, 360],
[2, 4, 0.05811, 0.17632, 0.034, 9900, 0, 0, 0, 0, 1, -360, 360],
[2, 5, 0.05695, 0.17388, 0.0346, 9900, 0, 0, 0, 0, 1, -360, 360],
[3, 4, 0.06701, 0.17103, 0.0128, 9900, 0, 0, 0, 0, 1, -360, 360],
[4, 5, 0.01335, 0.04211, 0, 9900, 0, 0, 0, 0, 1, -360, 360],
[4, 7, 0, 0.20912, 0, 9900, 0, 0, 0.978, 0, 1, -360, 360],
[4, 9, 0, 0.55618, 0, 9900, 0, 0, 0.969, 0, 1, -360, 360],
[5, 6, 0, 0.25202, 0, 9900, 0, 0, 0.932, 0, 1, -360, 360],
[6, 11, 0.09498, 0.1989, 0, 9900, 0, 0, 0, 0, 1, -360, 360],
[6, 12, 0.12291, 0.25581, 0, 9900, 0, 0, 0, 0, 1, -360, 360],
[6, 13, 0.06615, 0.13027, 0, 9900, 0, 0, 0, 0, 1, -360, 360],
[7, 8, 0, 0.17615, 0, 9900, 0, 0, 0, 0, 1, -360, 360],
[7, 9, 0, 0.11001, 0, 9900, 0, 0, 0, 0, 1, -360, 360],
[9, 10, 0.03181, 0.0845, 0, 9900, 0, 0, 0, 0, 1, -360, 360],
[9, 14, 0.12711, 0.27038, 0, 9900, 0, 0, 0, 0, 1, -360, 360],
[10, 11, 0.08205, 0.19207, 0, 9900, 0, 0, 0, 0, 1, -360, 360],
[12, 13, 0.22092, 0.19988, 0, 9900, 0, 0, 0, 0, 1, -360, 360],
[13, 14, 0.17093, 0.34802, 0, 9900, 0, 0, 0, 0, 1, -360, 360]])
ppc['bus'] = np.array([
[1, 3, 0, 0, 0, 0, 1, 1.06, 0, 0, 1, 1.06, 0.94, 232.4],
[2, 2, 21.7, 12.7, 0, 0, 1, 1.045, -4.98, 0, 1, 1.06, 0.94, 40.],
[3, 2, 94.2, 19, 0, 0, 1, 1.01, -12.72, 0, 1, 1.06, 0.94, 0.],
[4, 1, 47.8, -3.9, 0, 0, 1, 1.019, -10.33, 0, 1, 1.06, 0.94, 0.],
[5, 1, 7.6, 1.6, 0, 0, 1, 1.02, -8.78, 0, 1, 1.06, 0.94, 0.],
[6, 2, 11.2, 7.5, 0, 0, 1, 1.07, -14.22, 0, 1, 1.06, 0.94, 0.],
[7, 1, 0, 0, 0, 0, 1, 1.062, -13.37, 0, 1, 1.06, 0.94, 0.],
[8, 2, 0, 0, 0, 0, 1, 1.09, -13.36, 0, 1, 1.06, 0.94, 0.],
[9, 1, 29.5, 16.6, 0, 19, 1, 1.056, -14.94, 0, 1, 1.06, 0.94, 0.],
[10, 1, 9, 5.8, 0, 0, 1, 1.051, -15.1, 0, 1, 1.06, 0.94, 0.],
[11, 1, 3.5, 1.8, 0, 0, 1, 1.057, -14.79, 0, 1, 1.06, 0.94, 0.],
[12, 1, 6.1, 1.6, 0, 0, 1, 1.055, -15.07, 0, 1, 1.06, 0.94, 0.],
[13, 1, 13.5, 5.8, 0, 0, 1, 1.05, -15.16, 0, 1, 1.06, 0.94, 0.],
[14, 1, 14.9, 5, 0, 0, 1, 1.036, -16.04, 0, 1, 1.06, 0.94, 0.]])
if(mode=='Pow_full'): #This If is for the real power vector P
for r in range(14):
mat[r,0].value = ppc['bus'][r][2] +ppc['bus'][r][-1]
elif(mode=='B_full'): #This is the susceptance matrix
for r in range(14):
for c in range(14):
mat[r,c].value = 0.
for r in range(ppc['branch'].shape[0]):
fom = int(ppc['branch'][r][0])-1 #the from bus
tom = int(ppc['branch'][r][1])-1 #the to bus
mat[fom,tom].value = 1./ppc['branch'][r][3]
mat[tom,fom].value = 1./ppc['branch'][r][3]
for j in range(14):
mat[j,j].value = sum(mat[j])
else:
pass
return mat
谢谢
UPDATE 1
在功能dc_mats(mat,mode)
中,这部分代码引起了麻烦:
for j in range(14):
mat[j,j].value = sum(mat[j])
sum
返回数据类型instance
。但是,即使我注释了这段代码,我在定义的m.arrays
的优化部分仍然存在问题。
您的应用程序存在很多问题,因此我创建了一个使用随机初始值和矩阵初始值的简单应用程序。您的应用程序是一个线性方程组,因此应快速可靠地求解。您可以希望将您的问题特定信息填写到以下示例中。优化器调整u
的值,以将w
驱动到所需的设定点目标wref
。
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
m = GEKKO()
n = 14
B = np.ones((n,n))
H = np.ones(n)
wref = 0.5
u = m.Array(m.MV,n,lb=0,ub=1)
w = m.Array(m.Var,n)
theta = m.Array(m.Var,n)
P = np.dot(B,theta)
m.Equations([theta[i].dt()==w[i]-wref for i in range(n)])
m.Equations([w[i].dt()==(u[i]-P[i])/(2*H[i]) for i in range(n)])
[m.Minimize((w[i]-wref)**2) for i in range(n)]
m.time = np.linspace(0,5)
for i in range(n):
u[i].STATUS = 1
w[i].value = np.random.rand()
theta[i].value = np.random.rand()
m.options.IMODE = 6
m.options.SOLVER = 1
m.solve()
fig, (ax1,ax2,ax3) = plt.subplots(3,1)
for i in range(n):
ax1.plot(m.time,u[i].value)
ax2.plot(m.time,w[i].value)
ax3.plot(m.time,theta[i].value)
ax1.set_ylabel('u')
ax2.set_ylabel('w')
ax3.set_ylabel('theta')
ax2.plot([0,max(m.time)],[wref,wref],'k--',lw=3,label='Target')
ax2.legend()
ax3.set_xlabel('time')
plt.show()
我建议您查看类似的tutorial applications (see number 17 on MPC)或Machine Learning and Dynamic Optimization course上的应用程序。感谢您分享这个有趣的应用程序。