Python Gekko step-wise variable in a time domain problem with differential equations

I am coding an MPC problem in Python Gekko to heat a building on a winter day (for now I am working with a simple building with 2 zones only). The building is expressed with a Resistor-Capacitor (RC) model and the objective is to minimize a combination of maximum power and the total energy used. The RC model is expressed through differential equations. It is already working correctly by using the temperature of the zones as manipulated variables (MV) to minimize the objective.

But to increase the accuracy and integrate it better in the building I have now implemented a PID inside the MPC problem (because the building is using PID thermostats and I am using MPC only to feed better setpoints to the PID controllers). This means that the MPC will cahnge the PID setpoints to minimize the objective function. In this case the MVs are not the zones' temperatures but the PID setpoints. The zones temperatures will try to follow the setpoints through the PID control.

This is working well so far but I need the PID setpoints to change only every 2 hours instead of 15 minutes (which is the timestep I am using). The MPC must use 15 minutes timesteps as 2 hours would be way too large. This means that the MV should have step-wise shape with values that can change only every 8 timesteps. I am trying to address this with m.Array variables of the Setpoints (SP) (length of 97 as the m.time) to try to impose 12 different Fixed Values (FV) every 8 timesteps (SP[0:8] = FV_1 , SP[8:16] = FV_2 etc.) but I am not managing to solve my issue. I searched within multiple posts here in Stack Overflow but I am not finding a similar case.

The working code with the PID in the MPC is:

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. <= df.index.hour) & (df.index.hour < 18) # comfort schedules
heating_SP_unoccupied = 19.0 # Setpoints when occupied/unoccupied
heating_SP_occupied = 21.0
df["heating_SP"] = heating_SP_unoccupied
df.loc[df["is_occupied"], "heating_SP"] = heating_SP_occupied # Update for occupied setpoint
heating_SP = df["heating_SP"].values
T_comfort_21 = m.Param(value=heating_SP)

# PID parameters
Kc = 30000/3600             # controller gain
tauI = 1000                 # controller reset time
tauD = 0                    # derivative constant
Q1_0 = m.Const(value=0.0)   # OP bias Zone 1
Q2_0 = m.Const(value=0.0)   # OP bias Zone 2

# change solver options
m.options.IMODE = 6 # MPC

# parameters from building training
for i in params_list:
    locals()[i] = m.Const(value = locals()[i+'_results'][-1])

# initial conditions
T_wall_int_initial = np.ones(len(dates_mpc))*20
initial_thermal_power = U_test[:,0][0] + U_test[:,1][0]

# state variables (temperature)
T_1 = m.SV(value=X_test[:,0][0]);
T_2 = m.SV(value=X_test[:,1][0]);
T_wall_int = m.SV(value=T_wall_int_typique);

# manipulated variables (PID thermostat setpoints)
SP_1 = m.MV(value=heating_SP[0],fixed_initial=False); SP_1.STATUS = 1   # set point
SP_2 = m.MV(value=heating_SP[0],fixed_initial=False); SP_2.STATUS = 1   # set point
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

# controlled variables (heating)
Q_1 = m.Var(value=U_test[:,0][0],lb=0);
Q_2= m.Var(value=U_test[:,1][0],lb=0);
thermal_power = m.Var(value=initial_thermal_power)
max_power = m.FV(value=initial_thermal_power); max_power.STATUS = 1

# uncontrolled variables (weather and internal heat gains)
T_ext = m.Param(value=W_test[:,0])
Tground_1 = m.Param(value=W_test[:,1])
Tground_2 = m.Param(value=W_test[:,2])
Q_sol_1 = m.Param(value=W_test[:,3])
Q_sol_2 = m.Param(value=W_test[:,4])
Q_sol_roof = m.Param(value=W_test[:,5])
Q_int_1 = m.Param(value=W_test[:,6])
Q_int_2 = m.Param(value=W_test[:,7])
wind_sq = m.Param(value=W_test[:,8])

# building equations (RC grey box)
m.Equation(C_1 *1e6* T_1.dt() == (Q_1 + lambda_1 * Q_int_1 + sigma_1 * Q_sol_1 + sigma_roof * Q_sol_roof + alpha_1 * wind_sq * (T_ext - T_1) \
                    + Uext_1*(T_ext - T_1) + Uroof_1*(T_ext - T_1) + Ug_1*(Tground_1 - T_1) \
                    + U_wall_int*(T_wall_int - T_1) + U_12*(T_2- T_1)  ))

m.Equation(C_2 *1e6* T_2.dt() == (Q_2 + lambda_2 * Q_int_2 + sigma_2 * Q_sol_2 + sigma_roof * Q_sol_roof + alpha_2 * wind_sq * (T_ext - T_2) \
                    + Uext_2*(T_ext - T_2) + Uroof_2*(T_ext - T_2) + Ug_2*(Tground_2 - T_2) \
                    + U_wall_int*(T_wall_int - T_2) + U_12*(T_1- T_2) ))

m.Equation(C_wall_int *1e6* T_wall_int.dt() == U_wall_int*(T_1 - T_wall_int) + U_wall_int*(T_2 - T_wall_int))

# PID thermostats
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)

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

What I was trying to do is (i report here only the parts i modified):

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

I first tried to define the PID Setpoints as only m.Var but then i could not use SP_1[i] as it was saying "invalid index to scalar variable" so that's why I am using an Array and then I am trying to inpose something like SP_1[0:9] = an FV parameter that can change then SP_1[9:17] = another FV parameter that can change and so on, so that SP_1 is a step-wise variable. But now the m.Intermediate line err_1 = m.Intermediate(SP_1-T_1) gives this error: @error: Intermediate Definition Error: Intermediate variable with no equality (=) expression . I guess that it's because now SP_1 needs to be indexed there (since it's an m.Array) but then I cannot do the same with T_1 because that one needs to stay as State Variable (SV) because it is subject to differential equations and then also the integral error of the PID is subject to a differential equation. So I am basically stuck and don't know what to do.

What I wanted to ask is: Is it possible to implement the step-wise variable I would like to implement? and how could I do it?

Thanks to everyone that will answer to this post. If the question is not clear I can try to clarify it. Also, if it would be useful to have the input data to compute the code I could provide it.

Try using an option for the MV type that specifies the number of time steps between each allowable movement of the Manipulated Variable.

SP_1.MV_STEP_HOR=8
SP_2.MV_STEP_HOR=8

There is more information on MPC tuning options in the Dynamic Optimization Course and in the documentation.

A global option can be set for all MVs (e.g. m.options.MV_STEP_HOR=8) or they can be set individually (e.g. SP_1.MV_STEP_HOR=6, SP_1.MV_STEP_HOR=10). By default, the individual MVs are set to 0 when they use the global option.