CV-объект gekko использует только скалярное значение массива для SP, SPHI и SPLO, поэтому необходима некоторая модификация, чтобы оптимизатор рассмотрел будущие изменения уставки. Простое приложение MPC показывает, как заданные значения используются в Gekko.
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
m = GEKKO()
m.time = np.linspace(0,20,41)
p = m.MV(value=0, lb=0, ub=100) # Declare MV
p.STATUS = 1 # allow optimizer to change
p.DCOST = 0.1 # smooth MV response
p.DMAX = 10.0 # max move each cycle
v = m.CV(value=0) # Declare CV
v.STATUS = 1 # add CV to the objective
m.options.CV_TYPE = 2 # squared error
v.SP = 40 # set point
v.TR_INIT = 1 # set point trajectory
v.TAU = 5 # time constant of trajectory
m.Equation(10*v.dt() == -v + 2*p)
m.options.IMODE = 6 # control
m.solve(disp=False)
# get additional solution information
import json
with open(m.path+'//results.json') as f:
results = json.load(f)
plt.figure()
plt.subplot(2,1,1)
plt.plot(m.time,p.value,'b-',label='MV Optimized')
plt.legend()
plt.ylabel('Input')
plt.subplot(2,1,2)
plt.plot(m.time,results['v1.tr'],'k-',label='Reference Trajectory')
plt.plot(m.time,v.value,'r--',label='CV Response')
plt.ylabel('Output')
plt.xlabel('Time')
plt.legend(loc='best')
plt.show()
Существует одно заданное значение, которое является целью 40 иэталонная траектория определяется с постоянной времени 5. MPC не может точно следовать эталонной траектории из-за ограничения скорости изменения с DMAX=10. Есть два варианта, если вы хотите, чтобы оптимизатор знал о будущих изменениях заданного значения.
Вариант 1. Не используйте CV, используйте параметр прямой связи
Если вывам не нужна эталонная траектория, и все в порядке с квадратом целевой ошибки, тогда самый простой способ спроецировать будущие изменения уставки - это определить вашу собственную цель MPC с вектором параметров прямой связи значений уставок. Пример проблемы показывает, что оптимизатор ожидает изменения заданного значения и активно перемещается до следующего изменения заданного значения, чтобы минимизировать общую сумму квадратов ошибки. Это может быть желательно во многих обстоятельствах, но может быть нежелательно при изготовлении, когда происходят изменения качества продукта, и окончание производственной кампании должно быть в спецификации до производства переходного материала.
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
m = GEKKO()
m.time = np.linspace(0,20,41)
p = m.MV(value=0, lb=0, ub=100) # Declare MV
p.STATUS = 1 # allow optimizer to change
p.DCOST = 0.1 # smooth MV response
p.DMAX = 10.0 # max move constraint
v = m.Var(value=0)
sp = np.ones(41)*40
sp[20:] = 60
s = m.Param(value=sp)
m.Obj((s-v)**2)
m.Equation(10*v.dt() == -v + 2*p)
m.options.IMODE = 6 # control
m.solve(disp=False)
plt.figure()
plt.subplot(2,1,1)
plt.plot(m.time,p.value,'b-',label='MV Optimized')
plt.legend()
plt.ylabel('Input')
plt.subplot(2,1,2)
plt.plot(m.time,sp,'k-',label='Setpoint')
plt.plot(m.time,v.value,'r--',label='CV Response')
plt.ylabel('Output')
plt.xlabel('Time')
plt.legend(loc='best')
plt.show()
Вариант 2: Ошибка в качестве CV
Если желательно использовать опорную траекторию и встроенные опции CV Gekko, тогда можно выбратьопределите новую переменную ошибки e и управляйте ею. Переменная ошибки всегда имеет заданное значение, равное нулю, а заданное значение прямой связи реализовано как параметр прямой связи.
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
m = GEKKO()
m.time = np.linspace(0,20,41)
p = m.MV(value=0, lb=0, ub=100) # Declare MV
p.STATUS = 1 # allow optimizer to change
p.DCOST = 0.1 # smooth MV response
p.DMAX = 10.0 # max move constraint
v = m.Var(value=0)
sp = np.ones(41)*40
sp[20:] = 60
s = m.Param(value=sp)
e = m.CV(value=0) # Declare CV
e.STATUS = 1 # add CV to the objective
m.options.CV_TYPE = 2 # squared error
e.SP = 0 # set point
e.TR_INIT = 1 # error trajectory
e.TAU = 5 # time constant of trajectory
m.Equation(e==s-v)
m.Equation(10*v.dt() == -v + 2*p)
m.options.IMODE = 6 # control
m.solve(disp=False)
plt.figure()
plt.subplot(2,1,1)
plt.plot(m.time,p.value,'b-',label='MV Optimized')
plt.legend()
plt.ylabel('Input')
plt.subplot(2,1,2)
plt.plot(m.time,sp,'k-',label='Setpoint')
plt.plot(m.time,v.value,'r--',label='CV Response')
plt.ylabel('Output')
plt.xlabel('Time')
plt.legend(loc='best')
plt.show()