使用 GEKKO 建立电力系统模型
Power system model setup with GEKKO
我正在尝试使用 GEKKO. In specific using MPC to the IEEE 14 bus 测试用例优化电力系统。
系统包含 14 条总线，模型由状态变量 theta 和 omega（分别为发电机的功率角和旋转频率）和代数方程组成。
代数方程对应于直流近似，即潮流方程的线性近似。在下面的等式中，w^ref 是 60 Hz 的目标电频率。 P 是每条总线上的有功功率向量。 B是系统的电纳矩阵。 u是对应于发电机机械功率的控制输入。
objective是通过操纵机械功率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。
你的应用程序有很多问题，所以我创建了一个更简单的应用程序，它使用随机初始值和矩阵。您的应用程序是一个线性方程组，因此它应该可以快速可靠地求解。您可以在下面的示例中填写您的问题特定信息。优化器调整 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) or the applications on the Machine Learning and Dynamic Optimization course。感谢分享这个有趣的应用程序。
我正在尝试使用 GEKKO. In specific using MPC to the IEEE 14 bus 测试用例优化电力系统。
系统包含 14 条总线，模型由状态变量 theta 和 omega（分别为发电机的功率角和旋转频率）和代数方程组成。
代数方程对应于直流近似，即潮流方程的线性近似。在下面的等式中，w^ref 是 60 Hz 的目标电频率。 P 是每条总线上的有功功率向量。 B是系统的电纳矩阵。 u是对应于发电机机械功率的控制输入。
objective是通过操纵机械功率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。
你的应用程序有很多问题，所以我创建了一个更简单的应用程序，它使用随机初始值和矩阵。您的应用程序是一个线性方程组，因此它应该可以快速可靠地求解。您可以在下面的示例中填写您的问题特定信息。优化器调整 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) or the applications on the Machine Learning and Dynamic Optimization course。感谢分享这个有趣的应用程序。