Power system model setup with GEKKO

Question

I'm trying to optimize a power system with GEKKO. In specific using MPC to the IEEE 14 bus test case.

The system comprises 14 buses and the model consists of the state variables theta and omega (power angle and spinning frequency of the generator respectively) and the algebraic eqs.

The algebraic equations correspond to the DC approximation, i.e., linear approximation of the power flow equation. In the eqs below w^ref is the target electrical frequency of 60 Hz. P is the vector of real power at each bus. B is the susceptance matrix of the system. u is the control input which corresponds to the mechanical power of the generators.

The objective is to reach w^ref and power angles close to 0 for each bus by manipulating the mechanical power u.

The error I'm getting is (the code is below):

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()

My question is where is the coding wrong?

I have two function dc_opf() and dc_mats(mat,mode). The former is where the optimization happens. The latter is a helper function to populate the P and B matrices.

My code is:

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

Thanks

UPDATE 1

In the function dc_mats(mat,mode) This part of the code is causing troubles:

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

The sum is returning a datatype instance. But even if I comment this piece of code I still have issues in the optimization part with the m.arrays I'm defining.

user_2366941 · Answer

There were quite a few issues with your application so I created a simpler application that uses random initial values and ones for the matrices. Your application is a linear system of equations so it should solve quickly and reliably. You can hopefully fill in your problem specific information into the example below. The optimizer adjusts the value of u to drive w to the desired setpoint target 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()

I recommend that you look at similar tutorial applications (see number 17 on MPC) or the applications on the Machine Learning and Dynamic Optimization course. Thanks for sharing this interesting application.

Power system model setup with GEKKO

1 Answers