Implementing a MPC for regulation of ventilation system based on energy costs

We are a group of three members who are writing a bachelor thesis together about implementing MPC in a shopping mall here in Tromsø. Our main goal is to optimize the indoor temperature in based on power consumption. The issue is declaring all of our input variables and declaring the object function correctly.

All the input variables and are pulled from a pyads server from our Bechoff IPC where we get our real-time sensor data. We are also pulling weather data where we get the outdoor temperature for the next 12 hours in an array.

Our equations for the temperatures in the ventilation system:

T5_max = T1 + n1*(T3-T1)
T6_max = T5_max + n2*((V*rho*C_water)*n2*(T8-T7))/V_flow
T2_max = T5_max + T6_max + (n3 * V*rho*C_air*(T5-T1))/(V_flow * rho * C)
T5_min = T6_min = T2_min = T1

n1, n2, and n3 are all constants. (efficiency rates). All the T values are real-time sensor data. V_flow is constant. Our equations for the temperature in the mall.

sum(T3)dT/dt = (T3 - (UA*(T3-T1))/V*rho*C_air + ((T2-T3) * V_flow/Volume) dt
T_room = T3 + (sum(T3)dT/dt)

UA and Volume are known constants. We are also getting the power prices for a 24-hour period and the amount of people in the building. These are not implemented in the equations yet since we just want to get a basic model running first.

So, based on these equations we want to set the optimal temperature for t5, t6, and t2 to get the T_room in a range of 18-24 degrees. In the most cost-efficient way.

Hope someone can help us with implementing this in an easy way.

#Væskebatteri
NomEff_E2 = m.Param(value=50)   

#Gjenvinner 
n1 = 0.8
#E1 = VolumFlow*rho_luft*Cv_luft*(T5-uteTemp)

# Parameters
Volum = m.Const(value=39000) #m^3
VolumFlow = m.Const(value=5) #m^3/s
rho_luft = m.Const(value=1.2) #kg/m^3
rho_vann = m.Const(value=1000) #kg/m^3
Cv_vann = m.Const(value=4.18) #kJ/kgK
Cv_luft = m.Const(value=1.006) #kJ/kgK
NomEff_El = m.Const(value=100) #KW
T1 = m.Const(value=-10)#C
UA_verdi = m.Const(value=100) #W/C

t5 = m.Var() 
t6 = m.Var() 
t3 = m.Var() 
t2 = m.Var()


t5.value = T1
t6.value = T1
t2.value = T1
t3.VALUE = 18

#lower 
t3.lb = 18
t5.lb = T1
t6.lb = T1
t2.lb = T1

#upper
t5.ub = T1 + n1*(t3-T1)
t6.ub = T1 + n1*(t3-T1) + NomEff_E2/(VolumFlow*rho_luft*Cv_luft)
t2.ub = T1 + n1*(t3-T1) + NomEff_E2/(VolumFlow*rho_luft*Cv_luft)
                        + NomEff_El/(VolumFlow*rho_luft*Cv_luft)
t3.ub = 25

m.Equation(t3.dt == (t3-(UA_verdi*(t3-T1)))/
    (Volum*rho_luft*Cv_luft)+((t2-t3)*(VolumFlow/Volum)))

E1 = m.Param(VolumFlow*rho_luft*Cv_luft*(t5-T1))
E2 = m.Param(VolumFlow*rho_luft*Cv_luft*(t6-t5))
E3 = m.Param(VolumFlow*rho_luft*Cv_luft*(t2-t6))

n2 = m.Const(value = 0.2)
n3 = m.Const(value = 1)
spottpris = m.Const(value = 0.7)
effektledd = m.Const(value=15)

m.Obj((-E1 + E2*n2 + E3*n3)*(spottpris+effektledd))

m.options.IMODE = 3
m.solve(disp=False)
print(t3.VALUE, t5.VALUE, t6.VALUE, t2.VALUE)

This is what we were working on yesterday, but we can clearly tell that this isn't the way. Here we tried to only simulate with constant values and the variable declaration may not be the same as the equations above. We tried to follow the APMonitor example.

Consider us stupid.

Edit: Added a figure for a better explanation of our problem

Edit 2: Updated equations and hopefully more describing.

The system operates sequentially using up available energy at each point before tapping of the next. (E1 -> E2 -> E3)
The objective is to find the optimal T2 temperature set point for the air handling unit that minimizes costs. The spotprice array will be an array with twelve constants.

I've spent some time in Tromsø, so this is cool work to see. Here are a few general guidelines that should help.

1. Bounds for GEKKO variables cannot include GEKKO equations. They must evaluate to specific values that don't change during the optimization.
2. If variable bounds must be dynamic or vary depending on other variables in the problem they must be formulated as equations rather than as bounds in the variable declarations. See reworked code example below.
3. You can use normal Python values for constants that are not time-dependent and are not being optimized. This often keeps the code a little cleaner and easier to work with.
4. Try and define variable upper and lower bounds when you declare the GEKKO varible. If they must be defined later use the UPPER and LOWER properties rather than ub and lb, which are not defined.
5. If the program hangs with no visible progress(what happened when I first ran your code), you may have an unreported error on the remote server (default solution method). Solving the problem locally (m = GEKKO(remote=False)) can help get around this. Docs
6. The correct syntax for time-derivatives is t3.dt() not t3.dt. This one trips up a lot of people and was a problem in your code. If you ever get a "Model Expression" error it's worth checking to make sure you get this correct.
7. If you run into "Solution Not Found", your problem is likely mathmatically infeasible. You can open the solution directory (m.open_folder()) if you are running with remote=False and look at the infeasibilities.txt to see which equations are infeasible. If you add "names" to your variables like I did in the code example below it can help with reading the infeasibilities file.

Below is a working example of the code you posted. I have changed some variables names to help those that may not be familiar with Norwegian as well as simplified some things. I also added a slack variable to your first equation and the objective function to make the problem feasible.

There's a lot to take in there, so certainly let me know in the comments if any of this does not make sense.

from gekko import GEKKO

m = GEKKO(remote=False)

# Normal Python values can be used if the value is not 
# being optimized and is not time-dependent
NomEff_E2 = 50
n1 = 0.8            # recycle fraction
n2 = 0.2
n3 = 1
V = 39000           # m^3
AirFlow = 5         # m^3/s
rho_air = 1.2       # kg/m^3
rho_water = 1000    # kg/m^3
Cv_water = 4.18     # kJ/kgK
Cv_air = 1.006      # kJ/kgK
NomEff_El = 100     # KW
T1 = -10            # C, UteTemp
UA = 100            # W/C

# These will probably become m.parameters when you use real data
spotPrice = m.Const(value = 0.7)
effektledd = m.Const(value=15)

# The bounds on variables must be fixed values (not equations)
t5 = m.Var(T1, lb=T1, name='t5') 
t6 = m.Var(T1, lb=T1, name='t6') 
t3 = m.Var(18, lb=18, ub=25, name='t3') 
t2 = m.Var(T1, lb=T1, name='t2')

slack = m.Var(fixed_initial=False)

m.Equation(t3.dt() == (t3-(UA*(t3-T1)))/(V*rho_air*Cv_air)+(t2-t3)*(AirFlow/V) + slack) 

# The upper bounds expressed as equations
m.Equation(t5 <= T1 + n1*(t3-T1))
m.Equation(t6 <= T1 + n1*(t3-T1) + NomEff_E2/(AirFlow*rho_air*Cv_air))
m.Equation(t2 <= T1 + n1*(t3-T1) + NomEff_E2/(AirFlow*rho_air*Cv_air)
                        + NomEff_El/(AirFlow*rho_air*Cv_air))

E1 = m.Param(AirFlow*rho_air*Cv_air*(t5-T1))
E2 = m.Param(AirFlow*rho_air*Cv_air*(t6-t5))
E3 = m.Param(AirFlow*rho_air*Cv_air*(t2-t6))

m.Obj((-E1 + E2*n2 + E3*n3)*(spotPrice+effektledd)) + slack**2

m.options.IMODE = 3
m.solve(disp=False)
print(t3.VALUE, t5.VALUE, t6.VALUE, t2.VALUE)