gekko - 使用 GEKKO 建立电力系统模型

Question

我正在尝试使用 GEKKO 优化电力系统.具体使用 MPC 到 IEEE 14 bus测试用例。

该系统包含 14 条总线，模型由状态变量 theta 和 omega(分别为发电机的功率角和旋转频率)和代数方程组成。

代数方程对应于直流近似，即功率流方程的线性近似。在下面的方程式中，w^ref 是 60 Hz 的目标电频率。 P 是每条总线上有功功率的矢量。 B 是系统的电纳矩阵。 u是发电机机械功率对应的控制输入。

目标是通过操纵机械功率 u 使每辆总线的 w^ref 和功率角接近 0。

我得到的错误是(代码如下):

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

谢谢

更新 1

在函数dc_mats(mat,mode)这部分代码引起了麻烦:

for j in range(14):
 mat[j,j].value = sum(mat[j])

sum 返回数据类型instance。但是，即使我评论了这段代码，我在定义的 m.arrays 的优化部分仍然存在问题。

Answer 1

您的应用程序存在很多问题，因此我创建了一个更简单的应用程序，该应用程序使用随机初始值和矩阵初始值。您的应用程序是一个线性方程组，因此它应该可以快速可靠地求解。您可以在下面的示例中填写您的问题特定信息。优化器调整 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 上的应用程序.感谢分享这个有趣的应用程序。