Пошаговая переменная Python Gekko в задаче во временной области с дифференциальными уравнениями ⇐ Python

[Anonymous]

Я пишу задачу MPC в Python Gekko для обогрева здания в зимний день (сейчас я работаю с простым зданием только с двумя зонами). Здание представлено моделью резистора-конденсатора (RC), и цель состоит в том, чтобы минимизировать сочетание максимальной мощности и общей используемой энергии. Модель RC выражается через дифференциальные уравнения. Он уже работает правильно, используя температуру зон в качестве управляемых переменных (MV) для минимизации цели.
Но чтобы повысить точность и лучше интегрировать ее в здание, я сейчас реализовал ПИД внутри проблемы MPC (поскольку в здании используются ПИД-термостаты, а я использую MPC только для подачи лучших уставок на ПИД-регуляторы). Это означает, что MPC изменит уставки ПИД-регулятора, чтобы минимизировать целевую функцию. В этом случае значениями MV являются не температуры зон, а уставки ПИД. Температуры в зонах будут пытаться следовать заданным значениям посредством ПИД-регулирования.
Пока это работает хорошо, но мне нужно, чтобы заданные значения ПИД менялись только каждые 2 часа, а не 15 минут (что временной шаг, который я использую). MPC должен использовать временные интервалы в 15 минут, поскольку 2 часа были бы слишком большими. Это означает, что MV должен иметь ступенчатую форму со значениями, которые могут меняться только каждые 8 ​​тактов. Я пытаюсь решить эту проблему с помощью переменных m.Array уставок (SP) (длина 97 в качестве m.time), чтобы попытаться ввести 12 различных фиксированных значений (FV) каждые 8 ​​временных шагов (SP[0:8] = FV_1 , SP[8:16] = FV_2 и т. д.), но мне не удается решить свою проблему. Я искал несколько сообщений здесь, в Stack Overflow, но не нашел подобного случая.
Рабочий код с PID в MPC:

Код: Выделить всё

m = GEKKO(remote=False) # create GEKKO model
n_days_mpc = 1
starting_datetime_mpc = pd.Timestamp("2015-02-02 00:00", tz="US/Eastern")
finishing_datetime_mpc = starting_datetime_mpc +  timedelta(hours = (24*n_days_mpc))
dates_mpc = pd.date_range(starting_datetime_mpc, finishing_datetime_mpc,freq='15min',tz="US/Eastern")
time_array = np.arange(0,len(dates_mpc)*timestep,timestep)
m.time = time_array

# import input data
df, location, solar_position, Q_solar_df = weather_data (starting_datetime_mpc, finishing_datetime_mpc)
U_test = Power_delivered[starting_datetime_mpc:finishing_datetime_mpc].values
X_test = Temp[starting_datetime_mpc:finishing_datetime_mpc].values
W_cols = ["T_ext", "Tground_1", "Tground_2",
"Q_sol_1", "Q_sol_2", "Q_sol_roof",
"Q_int_1", "Q_int_2", "wind_sq"]  # exogenous_columns
W_test = df[W_cols].values

# reference performances from TRNSYS
power_trnsys = U_test[:,0] + U_test[:,1]
Q_int_trnsys = W_test[:,7] + W_test[:,8]
Q_sol_direct_trnsys = W_test[:,3] + W_test[:,4]

# comfort schedule
df["is_occupied"] = (7.  = thermal_power)

m.Equation(SP_1 >= T_comfort_21)
m.Equation(SP_2 >= T_comfort_21)

m.Minimize(14.58*max_power + 5.03/100*thermal_power/ (60/minutes_per_timestep))
m.solve()

Я пытался сделать следующее (здесь я сообщаю только о тех частях, которые я изменил):

Код: Выделить всё

# manipulated variables (PID thermostat setpoints)
SP_1 = m.Array(m.Var, len(m.time))
SP_2 = m.Array(m.Var, len(m.time))
SP_1_values = m.Array(m.FV,12)
SP_2_values = m.Array(m.FV,12)
for SP in SP_1_values:
SP.STATUS = 1
for SP in SP_2_values:
SP.STATUS = 1
Intgl_1 = m.Var((heating_SP[0]-X_test[:,0][0])*tauI*Kc)    # integral of the error Zone 1
Intgl_2 = m.Var((heating_SP[0]-X_test[:,1][0])*tauI*Kc)    # integral of the error Zone 2

# PID thermostats
for i in (0,1,2,3,4,5,6,7,8):
m.Equation(SP_1[i] == SP_1_values[0])
m.Equation(SP_2[i] == SP_2_values[0])
for j in range(1, 12):
start_index = 8 * j + 1
indices = tuple(range(start_index, start_index + 8))
for i in indices:
m.Equation(SP_1[i] == SP_1_values[j])
m.Equation(SP_2[i] == SP_2_values[j])
err_1 = m.Intermediate(SP_1-T_1) # set point error
err_2 = m.Intermediate(SP_2-T_2) # set point error
m.Equation(Intgl_1.dt()==err_1) # integral of the error
m.Equation(Q_1 == Q1_0 + Kc*err_1 + (Kc/tauI)*Intgl_1 - (Kc*tauD)*T_1.dt())
m.Equation(Intgl_2.dt()==err_2) # integral of the error
m.Equation(Q_2 == Q2_0 + Kc*err_2 + (Kc/tauI)*Intgl_2 - (Kc*tauD)*T_2.dt())

m.Equation(thermal_power == Q_1+Q_2)
m.Equation(max_power >= thermal_power)

for i in range(len(m.time)):
m.Equation(SP_1[i] >= T_comfort_21)
m.Equation(SP_2[i] >= T_comfort_21)

Сначала я попытался определить уставки ПИД только как m.Var, но затем не смог использовать SP_1, поскольку он говорил «неверный индекс скалярной переменной», вот почему я используя массив, а затем я пытаюсь ввести что-то вроде SP_1[0:9] = параметр FV, который может измениться, затем SP_1[9:17] = другой параметр FV, который может измениться, и так далее, так что SP_1 является шагом мудрая переменная. Но теперь строка m.Intermediate err_1 = m.Intermediate(SP_1-T_1) выдает следующую ошибку: @error: Ошибка определения промежуточного значения: Промежуточная переменная без выражения равенства (=) . Я предполагаю, что это потому, что теперь SP_1 должен быть проиндексирован там (поскольку это m.Array), но тогда я не могу сделать то же самое с T_1, потому что он должен оставаться переменной состояния (SV), потому что он подчиняется дифференциальным уравнениям, а затем также интегральная ошибка ПИД определяется дифференциальным уравнением. Так что я практически застрял и не знаю, что делать.
Я хотел спросить: возможно ли реализовать пошаговую переменную, которую я хотел бы реализовать? и как я мог это сделать?
Спасибо всем, кто ответит на этот пост. Если вопрос не ясен, я могу попытаться его уточнить. Кроме того, если бы было полезно иметь входные данные для вычисления кода, я мог бы их предоставить.

Подробнее здесь: https://stackoverflow.com/questions/791 ... tial-equat