利用GEKKO建立电力系统模型

Question

我正试图用盖科优化电力系统。具体而言，将MPC用于IEEE 14总线测试用例。

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

代数方程对应于直流近似，即潮流方程的线性近似。在低于w^ref的eqs中，目标电频率为60 Hz。P是每条总线上真实功率的矢量。B是系统的电纳矩阵。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()

我建议您查看类似的教程应用程序(参见MPC上的第17号)或机器学习与动态优化课程上的应用程序。谢谢你分享这个有趣的应用程序。