Почему GEKKO не предоставляет оптимальные команды, хотя выходные данные не соответствуют эталонным?
К этому вопросу относится следующее: модель прогнозирующего управления с использованием GEKKO
Я пытаюсь применить MPC для поддержания температуры в помещении в определенном диапазоне, но GEKKO выдает мне нулевые команды, даже если выходные данные расходятся. Я запускаю исправленный код из моего предыдущего вопроса:
# Import library
import numpy as np
import pandas as pd
import time
from gekko import GEKKO
from numpy import array
K = array([[ 0.93705481, -12.24012156]])
p = {'a': array([[ 1.08945247],
[-0.00242145],
[-0.00245978],
[-0.00272713],
[-0.00295845],
[-0.00319119],
[-0.00343511],
[-0.00366243],
[-0.00394247],
[-0.06665054]]),
'b': array([[[-0.05160235, -0.01039767],
[ 0.00173511, -0.01552485],
[ 0.00174602, -0.01179519],
[ 0.00180031, -0.01052658],
[ 0.00186416, -0.00822121],
[ 0.00193947, -0.00570905],
[ 0.00202877, -0.00344507],
[ 0.00211395, -0.00146947],
[ 0.00223514, 0.00021945],
[ 0.03800987, 0.04243736]]]),
'c': array([0.0265903])}
# i have used only 200 mes of T_externel
T_externel = np.linspace(9.51,9.78,200)
m = GEKKO(remote=False)
m.y = m.Array(m.CV,1)
m.u = m.Array(m.MV,2)
m.arx(p,m.y,m.u)
# rename CVs
m.T = m.y[0]
# rename MVs
m.beta = m.u[1]
# distrubance
m.d = m.u[0]
# distrubance and parametres
m.d = m.Param(T_externel[0])
# lower,heigh bound for MV
TL = m.Param(value = 16)
TH = m.Param(value = 18)
# steady state initialization
m.options.IMODE = 1
m.solve(disp=False)
# set up MPC
m.d.value = T_externel
m.options.IMODE = 6 # MPC
m.options.CV_TYPE = 2 # the objective is an l2-norm (squared error)
m.options.NODES = 2 # Collocation nodes
m.options.SOLVER = 1 # APOPT
m.time = np.arange(0,len(T_externel)*300,300) # step time = 300s
# Manipulated variables
m.beta.STATUS = 1 # calculated by the optimizer
m.beta.FSTATUS = 1 # use measured value
m.beta.DMAX = 1.0 # Delta MV maximum step per horizon interval
m.beta.DCOST = 0.0 # Delta cost penalty for MV movement
m.beta.UPPER = 1.0 # Lower bound
m.beta.LOWER = 0.0
m.beta.MEAS = 0 # set u=0
# Controlled variables
m.T.STATUS = 1 # drive to set point
m.T.FSTATUS = 1 # receive measurement
m.T.SP = 17 # set point
TL.value = np.ones(len(T_externel))*16
TH.value = np.ones(len(T_externel))*18
m.T.value = 17 # Temprature starts at 17
for i in range(len(T_externel)):
m.solve(disp = False)
if m.options.APPSTATUS == 1:
# Retrieve new values
beta = m.beta.NEWVAL
else:
# Solution failed
beta = 0.0
import matplotlib.pyplot as plt
# Plot the results
plt.figure(figsize=(12,6))
plt.subplot(2,1,1)
plt.plot(m.time,m.T.value,'r-',label=r'$T_{int}$')
plt.plot(m.time,TL.value,'k--',label='Lower Bound')
plt.plot(m.time,TH.value,'k--',label='Upper Bound')
plt.ylabel('Temperature (°C)')
plt.legend()
plt.subplot(2,1,2)
plt.plot(m.time,m.beta.value,'b--',label=r'$\beta$')
plt.ylabel('optimal control')
plt.xlabel('Time (sec)')
plt.legend()
plt.show()
И есть ли выход и оптимальное управление, полученные GEKKO:
Мне нужно поддерживать температуру в комнате в определенном диапазоне.
Можно быстро попробовать уменьшить подавление движений: m.beta.DCOST = 0.0. Можете ли вы сделать код автономным, включив в скрипт значения p из m.sysid()? Вы можете закомментировать шаги идентификации системы, и код сможет работать без набора данных. Также T_external = np.linspace(9.51,9.78,200) должно быть адекватной заменой беспокойства. Убедитесь, что код выполняется без каких-либо дополнительных внешних значений.
—
31.07.2024 02:32
@ Я внес рекомендованные вами изменения, но результаты все равно остались прежними. Я вставил код с внесенными изменениями (я включил значения p и k и установил m.beta.DCOST = 0.0).
—
31.07.2024 10:03
Почему в Python есть оператор "pass"?
Оператор pass в Python - это простая концепция, которую могут быстро освоить даже новички без опыта программирования.
Некоторые методы, о которых вы не знали, что они существуют в Python
Python - самый известный и самый простой в изучении язык в наши дни. Имея широкий спектр применения в области машинного обучения, Data Science,...
Основы Python Часть I
Вы когда-нибудь задумывались, почему в программах на Python вы видите приведенный ниже код?
LeetCode - 1579. Удаление максимального числа ребер для сохранения полной проходимости графа
Алиса и Боб имеют неориентированный граф из n узлов и трех типов ребер:
Оптимизация кода с помощью тернарного оператора Python
И последнее, что мы хотели бы показать вам, прежде чем двигаться дальше, это
Советы по эффективной веб-разработке с помощью Python
Как веб-разработчик, Python может стать мощным инструментом для создания эффективных и масштабируемых веб-приложений.
1
3
50
1
Перейти к ответу
Данный вопрос помечен как решенный
Ответы 1
Ответ принят как подходящий
Прирост обозначен как K = array([[ 0.93705481, -12.24012156]]), поэтому увеличение beta на +1 приводит к уменьшению T на -12,24. Увеличение T_ext на +1 приводит к увеличению T на +0,937. Стационарное значение T составляет 13,32, поэтому для получения начального значения 17 обычной практикой является создание аддитивного (или мультипликативного значения смещения, которое корректирует начальное значение. Добавляется дополнительная переменная Tb, хотя Gekko делает это внутренне.
При исследовании реакции контроллера хорошим первым шагом будет подтверждение выигрыша модели между MV и CV с помощью пошагового отклика. Вот ответ на шаг beta.
# Import library
import numpy as np
import pandas as pd
import time
from gekko import GEKKO
from numpy import array
K = array([[ 0.93705481, -12.24012156]])
p = {'a': array([[ 1.08945247],
[-0.00242145],
[-0.00245978],
[-0.00272713],
[-0.00295845],
[-0.00319119],
[-0.00343511],
[-0.00366243],
[-0.00394247],
[-0.06665054]]),
'b': array([[[-0.05160235, -0.01039767],
[ 0.00173511, -0.01552485],
[ 0.00174602, -0.01179519],
[ 0.00180031, -0.01052658],
[ 0.00186416, -0.00822121],
[ 0.00193947, -0.00570905],
[ 0.00202877, -0.00344507],
[ 0.00211395, -0.00146947],
[ 0.00223514, 0.00021945],
[ 0.03800987, 0.04243736]]]),
'c': array([0.0265903])}
# i have used only 200 mes of T_externel
T_externel = np.linspace(9.51,9.78,200)
m = GEKKO(remote=True)
m.y = m.Array(m.CV,1)
m.u = m.Array(m.MV,2)
m.arx(p,m.y,m.u)
# rename CVs
m.T = m.y[0]
# rename MVs
m.beta = m.u[1]
# distrubance
m.d = m.u[0]
# distrubance and parametres
m.d = m.Param(T_externel[0])
m.bias = m.Param(0)
m.Tb = m.CV()
m.Equation(m.Tb==m.T+m.bias)
# steady state initialization
m.options.IMODE = 1
m.solve(disp=True)
# set up MPC
#m.d.value = T_externel
m.options.IMODE = 6 # MPC
m.options.CV_TYPE = 1 # the objective is an l1-norm (region)
m.options.NODES = 2 # Collocation nodes
m.options.SOLVER = 3 # IPOPT
m.time = np.arange(0,len(T_externel)*300,300) # step time = 300s
# Manipulated variables
m.beta.STATUS = 0 # calculated by the optimizer
m.beta.FSTATUS = 0 # use measured value
m.beta.DMAX = 0.2 # Delta MV maximum step per horizon interval
m.beta.DCOST = 0.0 # Delta cost penalty for MV movement
m.beta.UPPER = 1.0 # Upper bound
m.beta.LOWER = 0.0 # Lower bound
m.beta.MEAS = 0.0 # Measured value
# step test
m.beta.value = np.zeros_like(m.time)
m.beta.value[20:] = 1
# Controlled variables
m.Tb.STATUS = 1 # drive to set point
m.Tb.FSTATUS = 1 # receive measurement
m.Tb.SPHI = 21 # set point high level
m.Tb.SPLO = 19 # set point low level
T_MEAS = 17 # Temperature starts at 17
m.Tb.value = T_MEAS
m.bias.value = T_MEAS - m.T.value[0]
m.solve(disp=True)
if m.options.APPSTATUS == 1:
# Retrieve new values
beta = m.beta.NEWVAL
else:
# Solution failed
beta = 0.0
import matplotlib.pyplot as plt
# Plot the results
plt.figure(figsize=(8,3.5))
plt.subplot(2,1,1)
plt.plot(m.time,m.Tb.value,'r-',label=r'$T_{biased}$')
plt.plot(m.time,m.T.value,'r--',label=r'$T_{unbiased}$')
plt.plot(m.time,m.d.value,'b:',label=r'T_{external}')
plt.plot([0,m.time[-1]],[m.Tb.SPHI,m.Tb.SPHI],'k--',label='Upper Bound')
plt.plot([0,m.time[-1]],[m.Tb.SPLO,m.Tb.SPLO],'k--',label='Lower Bound')
plt.ylabel('Temperature (°C)')
plt.legend(); plt.grid()
plt.subplot(2,1,2)
plt.step(m.time,m.beta.value,'b--',label=r'$\beta$')
plt.ylabel('optimal control')
plt.xlabel('Time (sec)')
plt.legend(); plt.grid()
plt.savefig('results.png',dpi=300)
plt.show()
Следующим шагом для исследования реакции контроллера является добавление возмущения и установка соответствующего диапазона уставки. Контроллер переключается на m.options.CV_TYPE=1, чтобы использовать цель l1-norm, которая дает опции m.Tb.SPHI и m.Tb.SPLO для указания диапазона.
# Import library
import numpy as np
import pandas as pd
import time
from gekko import GEKKO
from numpy import array
K = array([[ 0.93705481, -12.24012156]])
p = {'a': array([[ 1.08945247],
[-0.00242145],
[-0.00245978],
[-0.00272713],
[-0.00295845],
[-0.00319119],
[-0.00343511],
[-0.00366243],
[-0.00394247],
[-0.06665054]]),
'b': array([[[-0.05160235, -0.01039767],
[ 0.00173511, -0.01552485],
[ 0.00174602, -0.01179519],
[ 0.00180031, -0.01052658],
[ 0.00186416, -0.00822121],
[ 0.00193947, -0.00570905],
[ 0.00202877, -0.00344507],
[ 0.00211395, -0.00146947],
[ 0.00223514, 0.00021945],
[ 0.03800987, 0.04243736]]]),
'c': array([0.0265903])}
# i have used only 200 mes of T_externel
T_externel = np.linspace(9.51,9.78,200)
m = GEKKO(remote=True)
m.y = m.Array(m.CV,1)
m.u = m.Array(m.MV,2)
m.arx(p,m.y,m.u)
# rename CVs
m.T = m.y[0]
# rename MVs
m.beta = m.u[1]
# distrubance
m.d = m.u[0]
# distrubance and parametres
m.d = m.Param(T_externel[0])
m.bias = m.Param(0)
m.Tb = m.CV()
m.Equation(m.Tb==m.T+m.bias)
# steady state initialization
m.options.IMODE = 1
m.solve(disp=True)
# set up MPC
m.d.value = T_externel
m.options.IMODE = 6 # MPC
m.options.CV_TYPE = 1 # the objective is an l1-norm (region)
m.options.NODES = 2 # Collocation nodes
m.options.SOLVER = 3 # IPOPT
m.time = np.arange(0,len(T_externel)*300,300) # step time = 300s
# Manipulated variables
m.beta.STATUS = 1 # calculated by the optimizer
m.beta.FSTATUS = 0 # use measured value
m.beta.DMAX = 1.0 # Delta MV maximum step per horizon interval
m.beta.DCOST = 0.0 # Delta cost penalty for MV movement
m.beta.UPPER = 1.0 # Upper bound
m.beta.LOWER = 0.0 # Lower bound
# Controlled variables
m.Tb.STATUS = 1 # drive to set point
m.Tb.FSTATUS = 0 # receive measurement
m.Tb.SPHI = 15 # set point high level
m.Tb.SPLO = 13 # set point low level
m.Tb.WSPHI = 100 # set point high priority
m.Tb.WSPLO = 100 # set point low priority
T_MEAS = 17 # Temperature starts at 17
m.Tb.value = T_MEAS
m.bias.value = T_MEAS - m.T.value[0]
m.solve(disp=True)
if m.options.APPSTATUS == 1:
# Retrieve new values
beta = m.beta.NEWVAL
else:
# Solution failed
beta = 0.0
import matplotlib.pyplot as plt
# Plot the results
plt.figure(figsize=(8,3.5))
plt.subplot(2,1,1)
plt.plot(m.time,m.Tb.value,'r-',label=r'$T_{biased}$')
plt.plot(m.time,m.T.value,'r--',label=r'$T_{unbiased}$')
plt.plot(m.time,m.d.value,'g:',label=r'$T_{ext}$')
plt.plot([0,m.time[-1]],[m.Tb.SPHI,m.Tb.SPHI],'k--',label='Upper Bound')
plt.plot([0,m.time[-1]],[m.Tb.SPLO,m.Tb.SPLO],'k--',label='Lower Bound')
plt.ylabel('Temperature (°C)')
plt.legend(loc=1); plt.grid()
plt.subplot(2,1,2)
plt.step(m.time,m.beta.value,'b--',label=r'$\beta$')
plt.ylabel('optimal control')
plt.xlabel('Time (sec)')
plt.legend(loc=1); plt.grid()
plt.savefig('results.png',dpi=300)
plt.show()
Контроллер удерживает значение Tb в пределах заданного диапазона, регулируя beta.
Также возможно сделать значение beta целым числом, если система охлаждения имеет только состояние включения/выключения вместо непрерывных значений от 0 до 1.
# Import library
import numpy as np
import pandas as pd
import time
from gekko import GEKKO
from numpy import array
K = array([[ 0.93705481, -12.24012156]])
p = {'a': array([[ 1.08945247],
[-0.00242145],
[-0.00245978],
[-0.00272713],
[-0.00295845],
[-0.00319119],
[-0.00343511],
[-0.00366243],
[-0.00394247],
[-0.06665054]]),
'b': array([[[-0.05160235, -0.01039767],
[ 0.00173511, -0.01552485],
[ 0.00174602, -0.01179519],
[ 0.00180031, -0.01052658],
[ 0.00186416, -0.00822121],
[ 0.00193947, -0.00570905],
[ 0.00202877, -0.00344507],
[ 0.00211395, -0.00146947],
[ 0.00223514, 0.00021945],
[ 0.03800987, 0.04243736]]]),
'c': array([0.0265903])}
# i have used only 200 mes of T_externel
T_externel = np.linspace(9.51,9.78,200)
m = GEKKO(remote=True)
m.y = m.Array(m.CV,1)
m.u = m.Array(m.MV,2)
m.arx(p,m.y,m.u)
# rename CVs
m.T = m.y[0]
# rename MVs
m.beta = m.u[1]
m.free(m.beta)
m.bint = m.MV(0,lb=0,ub=1,integer=True)
m.Equation(m.beta==m.bint)
# distrubance
m.d = m.u[0]
# distrubance and parametres
m.d = m.Param(T_externel[0])
m.bias = m.Param(0)
m.Tb = m.CV()
m.Equation(m.Tb==m.T+m.bias)
# steady state initialization
m.options.IMODE = 1
m.solve(disp=True)
# set up MPC
m.d.value = T_externel
m.options.IMODE = 6 # MPC
m.options.CV_TYPE = 1 # the objective is an l1-norm (region)
m.options.NODES = 2 # Collocation nodes
m.options.SOLVER = 3 # IPOPT
m.time = np.arange(0,len(T_externel)*300,300) # step time = 300s
# Manipulated variables
m.bint.STATUS = 1 # calculated by the optimizer
m.bint.FSTATUS = 0 # use measured value
m.bint.DCOST = 0.0 # Delta cost penalty for MV movement
m.bint.UPPER = 1.0 # Upper bound
m.bint.LOWER = 0.0 # Lower bound
m.bint.MV_STEP_HOR = 10
m.bint.value = 0
# Controlled variables
m.Tb.STATUS = 1 # drive to set point
m.Tb.FSTATUS = 0 # receive measurement
m.Tb.SPHI = 15 # set point high level
m.Tb.SPLO = 13 # set point low level
m.Tb.WSPHI = 100 # set point high priority
m.Tb.WSPLO = 100 # set point low priority
T_MEAS = 17 # Temperature starts at 17
m.Tb.value = T_MEAS
m.bias.value = T_MEAS - m.T.value[0]
m.options.SOLVER = 1
m.solve(disp=True)
if m.options.APPSTATUS == 1:
# Retrieve new values
beta = m.beta.NEWVAL
else:
# Solution failed
beta = 0.0
import matplotlib.pyplot as plt
# Plot the results
plt.figure(figsize=(8,3.5))
plt.subplot(2,1,1)
plt.plot(m.time,m.Tb.value,'r-',label=r'$T_{biased}$')
plt.plot(m.time,m.T.value,'r--',label=r'$T_{unbiased}$')
plt.plot(m.time,m.d.value,'g:',label=r'$T_{ext}$')
plt.plot([0,m.time[-1]],[m.Tb.SPHI,m.Tb.SPHI],'k--',label='Upper Bound')
plt.plot([0,m.time[-1]],[m.Tb.SPLO,m.Tb.SPLO],'k--',label='Lower Bound')
plt.ylabel('Temperature (°C)')
plt.legend(loc=1); plt.grid()
plt.subplot(2,1,2)
plt.step(m.time,m.bint.value,'b--',label=r'$\beta_{int}$')
plt.ylabel('optimal control')
plt.xlabel('Time (sec)')
plt.legend(loc=1); plt.grid()
plt.savefig('results.png',dpi=300)
plt.show()
Время решения для двоичных переменных значительно больше, поэтому многие из них составляют приблизительно от 0,37 до 37% ВКЛ с некоторой постобработкой, например, выключение 3,7 минуты и включение 6,3 минуты в блоках с 10-минутными интервалами.
Iter: 71 I: 0 Tm: 16.99 NLPi: 10 Dpth: 4 Lvs: 87 Obj: 3.04E+02 Gap: 6.24E-01
Iter: 72 I: 0 Tm: 11.37 NLPi: 3 Dpth: 4 Lvs: 88 Obj: 3.04E+02 Gap: 6.24E-01
--Integer Solution: 3.04E+02 Lowest Leaf: 3.04E+02 Gap: 9.61E-09
Iter: 73 I: 0 Tm: 13.62 NLPi: 3 Dpth: 13 Lvs: 88 Obj: 3.04E+02 Gap: 9.61E-09
Successful solution
---------------------------------------------------
Solver : APOPT (v1.0)
Solution time : 1064.45290000000 sec
Objective : 304.259454634496
Successful solution
---------------------------------------------------
Большое спасибо, профессор @John Hedengern, за ваши объяснения и за то, что нашли время помочь мне решить мою проблему. Теперь я понимаю вещи гораздо лучше.
—
01.08.2024 18:27
Другие вопросы по теме
Проблема оптимизации расписания не знаю лучшего способа представить решение — Python/Gekko
Почему он предоставляет два разных вывода с помощью if2/if3?
Ограничения одновременного интервала и продолжительности с временными промежутками в Gekko
Ограничения по времени в Gekko
Работа с неоптимальными решениями от Gekko
Динамическая оптимизация с линейной интерполяцией по некоторым параметрам в GEKKO
Линейная интерполяция в ОДУ (ГЕККО)
Как изменить статус MV переменной GEKKO после решения устойчивого состояния?
Что вызывает эту синтаксическую ошибку Gekko: Ошибка в синтаксисе функциональной строки: отсутствует открывающая скобка?
Подгонка двух совокупностей к измерениям с использованием GEKKO, как оптимизировать первую точку данных
Похожие вопросы
Java sshtools сгенерировал подпись EDDSA, не совпадающую с сгенерированной подписью пикриптома Python
Как записать атрибуты каждого экземпляра в переменную класса, при этом родительский класс записывает экземпляры всех классов в одну и ту же переменную?
Создание динамического массива numpy с использованием двух существующих массивов
Функция C++ возвращает результат очень медленно, намного медленнее, чем функционально эквивалентный код Python
Добавление переменного отступа перед каждой строкой вывода журнала Python
Как извлечь регулярное выражение с переменной из строки в пандах?
Создайте PyObject с указателем на функцию C
Как аннотировать тип OrderedDict, который инициализируется литералами?
Как cplot() условная функция
Как получить DataFrame, разделить его на более мелкие сегменты и записать их в книгу Excel на разных листах
@ Джон, я забыл упомянуть, что на третьем рисунке показаны значения фиксированной наружной температуры.