Can someone help me get this Gekko model to work? It is about an optimization problem in which the best hourly control values for the compressor of a refrigeration system are to be found depending on a price vector, storage size and the cooling load for a 24-hour period.
I actually want to implement machine learning models into this optimization problem. However, so far it doesn't even work with test functions.
from gekko import GEKKO
from gekko import ML
from gekko.ML import Gekko_GPR, Gekko_SVR, Gekko_LinearRegression
# Initialize Gekko model
m = GEKKO(remote=False) # 'remote=False' for local execution
price_vec_list = [1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 2, 1, 1, 1] # Price vector
# Create a Gekko array of parameters
price_vec_gekko = m.Array(m.Param, len(price_vec_list))
for i, val in enumerate(price_vec_list):
price_vec_gekko[i].value = val
coolLoad = m.Const(3) # Cooling load
Q_Spc = m.Const(5) # Storage size
n = len(price_vec_gekko) # Directly use the length of the price vector
T_brine_in = m.Array(m.Param, n) # Brine inlet temperature
for i in range(n):
T_brine_in[i].VALUE = 25
ST_V = m.Array(m.MV, n) # Compressor control variable
for i in range(n):
ST_V[i].VALUE = 0.5
ST_V[i].LOWER = 0.43 # Correct property for lower limit
ST_V[i].UPPER = 1 # Correct property for upper limit
ST_V[i].status = 1 # Allow the solver to optimize this variable
coolLoad_sum_min = np.cumsum(np.full(n, coolLoad))
coolLoad_sum_max = coolLoad_sum_min + Q_Spc # Cold water storage cap
costs = m.Var(value=0) # Initialize cost variable as a Gekko variable
# Equations Object und Constraints
for i in range(n):
# Calculate the increment in costs directly within the equation
m.Equation(costs == costs + (2.5 * ST_V[i] + 0.01 * T_brine_in[i]) * price_vec_gekko[i])
# Define the Qo_constraint directly within the constraints
m.Equation(-6 * (ST_V[i] - 1) ** 2 + 6 + 0.01 * T_brine_in[i] < coolLoad_sum_max[i])
m.Equation(-6 * (ST_V[i] - 1) ** 2 + 6 + 0.01 * T_brine_in[i] > coolLoad_sum_min[i])
#for i in range(n):
#costs_increment = Gekko_LinearRegression(MLA_LR_P, m).predict([ST_V[i], T_brine_in[i]]) / 1000 * price_vec_gekko[i]
#Qo_constraint = Gekko_LinearRegression(MLA_LR_Q, m).predict([ST_V[i], T_brine_in[i]]) / 1000
## Update costs based on predictions and prices
#m.Equation(costs == costs + costs_increment)
## Constraints based on Qo_Gekko_Model, ensure these are properly defined as equations
#m.Equation(Qo_constraint < coolLoad_sum_max[i])
#m.Equation(Qo_constraint > coolLoad_sum_min[i])
m.Minimize(costs) # Minimize
m.options.SOLVER = 1 # APOPT Solver
m.options.IMODE = 3
# Solve model
m.solve(disp=True)
The output was:
APM Model Size
Each time step contains
Objects : 0
Constants : 2
Variables : 121
Intermediates: 0
Connections : 0
Equations : 73
Residuals : 73
Number of state variables: 73
Number of total equations: - 72
Number of slack variables: - 48
Degrees of freedom : -47
* Warning: DOF <= 0
Steady State Optimization with APOPT Solver
Iter Objective Convergence
0 1.93114E+05 9.03125E+00
1 -5.00483E+04 1.87500E-01
2 -6.39512E+11 2.43675E-01
.......
248 -2.14152E+14 2.43675E-01
249 -2.14919E+14 2.43675E-01
Iter Objective Convergence
250 -2.15686E+14 2.43675E-01
Maximum iterations
Solver : APOPT (v1.0)
Solution time : 2.1469 sec
Objective : -6.0710150204877865E+28
Unsuccessful with error code 0
Creating file: infeasibilities.txt
Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
and then I get an error at m.solve(disp=True):
@error: Solution Not Found
Exception Traceback (most recent call last)
Cell In[33], line 52
50 m.options.IMODE = 3
51 # Solve model
52 m.solve(disp=True)
54 # Output results
55 print(f'Optimized ST_V: {ST_V.VALUE[0]}')
File ~\anaconda3\lib\site-packages\gekko\gekko.py:2140, in GEKKO.solve(self, disp, debug, GUI, **kwargs)
2138 print("Error:", errs)
2139 if (debug >= 1) and record_error:
-> 2140 raise Exception(apm_error)
2142 else: #solve on APM server
2143 def send_if_exists(extension):
Exception: @error: Solution Not Found
I have already been able to solve the optimization problem with scipy. Maybe that helps to understand what the goal is. We came up with Gekko because Gekko should also be able to solve MINLP problems. We currently only modulate the compressor control value between the minimum value of 0.42 and the maximum value of 1. Switching off is not possible with scipy or with fmincon in matlab. In this code you can see how I solved it in scipy (different machine learning models can be loaded at the beginning. MLA_PR5_P and MLA_PR5_Q are polynomial ML models of the 5th degree):
import numpy as np
import pandas as pd
import joblib
from scipy.optimize import minimize
# Load models once
MLA_P_model = joblib.load('MLA_PR5_P.joblib')
MLA_Q_model = joblib.load('MLA_PR5_Q.joblib')
def predict_with_MLA_P(ST_V, T_brine_in, model):
# Use the already loaded model
input_data_df = pd.DataFrame({'ST_V': ST_V, 'T_brine_in': T_brine_in})
prediction = model.predict(input_data_df)
return prediction
def predict_with_MLA_Q(ST_V, T_brine_in, model):
# Use the already loaded model
input_data_df = pd.DataFrame({'ST_V': ST_V, 'T_brine_in': T_brine_in})
prediction = model.predict(input_data_df)
return prediction
def objective_function(ST_V, T_brine_in, price_vec):
Pel = predict_with_MLA_P(ST_V, T_brine_in, MLA_P_model)/1000
costs = np.sum(Pel * price_vec)
return costs
def constraint1(ST_V, T_brine_in, Qo_SOC_start, coolLoad, Q_Spc):
# Calculates the difference between Qo_sum and coolLoad_sum_max
Qo = predict_with_MLA_Q(ST_V, T_brine_in, MLA_Q_model)/1000
Qo_sum = Qo_SOC_start + np.cumsum(Qo)
coolLoad_sum_min = np.cumsum(np.full(len(ST_V), coolLoad))
coolLoad_sum_max = coolLoad_sum_min + Q_Spc
return coolLoad_sum_max - Qo_sum
def constraint2(ST_V, T_brine_in, Qo_SOC_start, coolLoad):
# Calculates the difference between coolLoad_sum_min and Qo_sum
Qo = predict_with_MLA_Q(ST_V, T_brine_in, MLA_Q_model)/1000
Qo_sum = Qo_SOC_start + np.cumsum(Qo)
coolLoad_sum_min = np.cumsum(np.full(len(ST_V), coolLoad))
return Qo_sum - coolLoad_sum_min
# Define parameters
n = 24
T_brine_in = np.random.randint(21, 40, size=n)
price_vec = [1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 2, 1, 1, 1]
Qo_SOC_start = 0
coolLoad = 5
Q_Spc = 5
ST_V_Start = np.ones(n) * 0.5 # Initial value, 1D vector with 24 rows
# Correction and definition of constraints as a list of dictionaries
constraints = [
{'type': 'ineq', 'fun': lambda ST_V: constraint1(ST_V, T_brine_in, Qo_SOC_start, coolLoad, Q_Spc)},
{'type': 'ineq', 'fun': lambda ST_V: constraint2(ST_V, T_brine_in, Qo_SOC_start, coolLoad)}
]
# Use in minimize function
result = minimize(
fun=objective_function,
x0=ST_V_Start,
args=(T_brine_in, price_vec),
constraints=constraints,
bounds=[(0.42, 1)] * n,
method='SLSQP'
)
print(result)
print(result.x)
message: Optimization terminated successfully
success: True
status: 0
fun: 64.15406185471795
x: [ 8.871e-01 8.892e-01 ... 7.418e-01 7.268e-01]
nit: 27
jac: [ 2.720e+00 2.622e+00 ... 2.631e+00 2.463e+00]
nfev: 694
njev: 27
[0.8871013 0.88922617 0.89714527 0.92713093 0.89713152 0.927064
0.89716403 0.92715987 0.90773437 0.7351564 0.70690274 0.69600369
0.71458654 0.69350165 0.6934931 0.74089652 0.54573353 0.543359
0.51353771 0.51842201 0.70690115 0.74303588 0.74178497 0.72679732]
Here are two plots for further understanding
As you can see, the compressor control value is high at a low prices (1) and low at a high prices (3). In the second plot you can see that the cold produced is within the limits between the cooling load and the cooling load + storage capacity. I also tried converting coolLoad_sum_min and coolLoad_sum_max to Gekko.arrays. That didn't help either.
It looks like you have an infeasibile constraint with the equation. You can look through the model structure and possible infeasible constraints by calling m.open_folder(), and comparing gk_modelX.apm and infeasibilities.txt. Printing out gekko variable names will also help identify which equations are infeasible.
This equation stood out to me as a problem:
The costs is just one variable, not an array, so this is trying to make the rest of the equation ((2.5 * ST_V[i] + 0.01 * T_brine_in[i]) * price_vec_gekko[i]) 0, which may not be feasible. Generally, use m.Equation for equality and inequality constraints; if you are defining a new variable/objective value (like costs), you can define it as an intermediate with m.Intermediate(). Also, by defining costs as m.Var(), you are defining it as a free variable that can be changed by the solver; I suggest initializing it as 0 and updating it like so;
After some tweaking, I was not able to get the code to solve. It looks like the lower limit is too restricting for how you calculate the Qo_Constraint; the max possible value with the parameters you provide is 6.25, while the minimum constraint rises well above that. Here is a slimmed down code that I used:
Hopefully this helps!