Question

Я пытаюсь оптимизировать систему питания с помощью GEKKO . В c используется MP C для теста IEEE 14 .

Система состоит из 14 шин, а модель состоит из переменных состояния theta и omega (угол мощности и частота вращения генератора соответственно) и алгебра c уравнения

Уравнения алгебры c соответствуют приближению D C, т. е. линейному приближению уравнения потока мощности. В уравнениях ниже w^ref - целевая электрическая частота 60 Гц. 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

Спасибо

ОБНОВЛЕНИЕ 1

В функции dc_mats(mat,mode) Эта часть кода вызывает проблемы:

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

sum возвращает тип данных instance. Но даже если я прокомментирую этот фрагмент кода, у меня все еще есть проблемы в части оптимизации с m.arrays, который я определяю.

John Hedengren · Answer 1 · 25 марта 2020

Было довольно много проблем с вашим приложением, поэтому я создал более простое приложение, которое использует случайные начальные значения и значения для матриц. Ваша заявка представляет собой линейную систему уравнений, поэтому она должна решаться быстро и надежно. Мы надеемся, что вы можете указать свою проблему, указав c информацию в следующем примере. Оптимизатор настраивает значение 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()

Я рекомендую вы смотрите на похожие учебные приложения (см. номер 17 на MP C) или приложения на Машинное обучение и Dynami c курс оптимизации . Спасибо за то, что поделились этим интересным приложением.

Настройка модели энергосистемы с GEKKO

Пожалуйста, войдите или зарегистрируйтесь чтобы ответить на этот вопрос.

1 Ответ

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Настройка модели энергосистемы с GEKKO

Пожалуйста, войдите или зарегистрируйтесь чтобы ответить на этот вопрос.

1 Ответ

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Похожие темы