I ran the following code in Excel and received I_g=[0,0,0], I_d=[14.8, 14.3, 10.3] and Objective=383.546, but I am having a problem to get a same results via GEKKO in Python. I used the four last constraints to enforce either I_g or I_d become 1 when their corresponding S_g or S_d is greater than zero, but both I_d and I_g are zero again. Maybe, my problem is the binary variables' definition in Gekko, because I expect the last four constraints check the following conditions: If S_g>0 then I_g=1 else I_g=0, If S_d>0 then I_d=1 else I_d=0. While S_d is greater than zero, but I_d is again 0. I appreciate your help.
from gekko import GEKKO
import pandas as pd
import sys
Inp_Data=pd.DataFrame([[10.8,10.6,15.0,0.2,0.0],[10.2,10.3,15.0,0.7,0.0],[10.1,10.4,15.0,4.9,0.0]])
col_name=list(Inp_Data)
row=len(Inp_Data)
P_D=Inp_Data.loc[:,col_name[0]].tolist()
P_I=Inp_Data.loc[:,col_name[1]].tolist()
Demand=Inp_Data.loc[:,col_name[2]].tolist()
Cap=300
FCV=[0.246,0.08145]
V_V=0.001
G_w=Inp_Data.loc[:,col_name[3]]
G_p=Inp_Data.loc[:,col_name[4]]
Big_M=100000000
epsilon=sys.float_info.epsilon
m = GEKKO(remote=False)
m.options.SOLVER=1
m.solver_options = ['minlp_gap_tol 1.0e-2',\
'minlp_maximum_iterations 10000',\
'minlp_max_iter_with_int_sol 500',\
'minlp_branch_method 1',\
'minlp_integer_leaves 2']
#Definition of the MIP variables
I_g=[m.Var(integer=True,lb=0,ub=1) for i in range(row)]
S_g=[m.Var(lb=0) for i in range(row)]
I_v=[m.Var(integer=True,lb=0,ub=1) for i in range(row)]
S_v=[m.Var(lb=0) for i in range(row)]
I_d=[m.Var(integer=True,lb=0,ub=1) for i in range(row)]
S_d=[m.Var(lb=0) for i in range(row)]
G_d=[m.Var(lb=0) for i in range(row)]
I_w=[m.Var(integer=True,lb=0,ub=1) for i in range(row)]
I_p=[m.Var(integer=True,lb=0,ub=1) for i in range(row)]
I_DA=[m.Var(integer=True,lb=0,ub=1)]
c = [None]*row
for t in range(row):
c[t] = (I_DA[0]*Demand[t]*P_D[t]+(1-I_DA[0])*S_v[t]*P_I[t])-(I_g[t]*P_I[t]*S_g[t]+I_d[t]*(FCV[0]*S_d[t]+FCV[1]*Cap)+V_V*S_v[t])
m.Equations(S_v[t]+S_g[t]+S_d[t]-(I_DA[0]*Demand[t]+(1-I_DA[0])*S_v[t])==0 for t in range(row))
m.Equations(S_v[t]==(G_w[t]+G_p[t]) for t in range(row))
m.Equations(I_d[t]+I_g[t]-I_DA[0]<=0 for t in range(row))
m.Equations(I_g[t]<(1+S_g[t]-epsilon) for t in range(row))
m.Equations(S_g[t]<(Big_M*I_g[t]) for t in range(row))
#
m.Equations(I_d[t]<(1+S_d[t]-epsilon) for t in range(row))
m.Equations(S_d[t]<(Big_M*I_d[t]) for t in range(row))
m.Maximize(m.sum(c))
m.solve(disp=True)
print('Results')
for t in range(row):
print('I_g{}={}, I_d{}={}'.format(t,str(I_g[t].value),t,str(I_d[t].value)))
print('Objective: ' + str(m.options.objfcnval))
Add 'minlp_integer_tol 1e-8' to the solver options to get the correct solution. The default is 'minlp_integer_tol 1e-2' meaning that a binary value can be 0.01 or 0.99 and is still considered an integer solution. With the new tolerance, a value of 1e-8 or 0.99999999 is considered an integer value. Here is the new solution with the tolerance applied:
I_g0=[0.0], I_d0=[1.0]
I_g1=[0.0], I_d1=[1.0]
I_g2=[0.0], I_d2=[1.0]
Objective: -383.546
Here is the complete program:
from gekko import GEKKO
import pandas as pd
import sys
Inp_Data=pd.DataFrame([[10.8,10.6,15.0,0.2,0.0],\
[10.2,10.3,15.0,0.7,0.0],\
[10.1,10.4,15.0,4.9,0.0]])
col_name=list(Inp_Data)
row=len(Inp_Data)
P_D=Inp_Data.loc[:,col_name[0]].tolist()
P_I=Inp_Data.loc[:,col_name[1]].tolist()
Demand=Inp_Data.loc[:,col_name[2]].tolist()
Cap=300
FCV=[0.246,0.08145]
V_V=0.001
G_w=Inp_Data.loc[:,col_name[3]]
G_p=Inp_Data.loc[:,col_name[4]]
Big_M=100000000
epsilon=sys.float_info.epsilon
m = GEKKO(remote=False)
m.options.SOLVER=1
m.solver_options = ['minlp_gap_tol 1.0e-2',\
'minlp_maximum_iterations 10000',\
'minlp_max_iter_with_int_sol 500',\
'minlp_branch_method 1',\
'minlp_integer_tol 1e-8', \
'minlp_integer_leaves 2']
#Definition of the MIP variables
I_g=[m.Var(integer=True,lb=0,ub=1) for i in range(row)]
S_g=[m.Var(lb=0) for i in range(row)]
I_v=[m.Var(integer=True,lb=0,ub=1) for i in range(row)]
S_v=[m.Var(lb=0) for i in range(row)]
I_d=[m.Var(integer=True,lb=0,ub=1) for i in range(row)]
S_d=[m.Var(lb=0) for i in range(row)]
G_d=[m.Var(lb=0) for i in range(row)]
I_w=[m.Var(integer=True,lb=0,ub=1) for i in range(row)]
I_p=[m.Var(integer=True,lb=0,ub=1) for i in range(row)]
I_DA=[m.Var(integer=True,lb=0,ub=1)]
c = [None]*row
for t in range(row):
c[t] = (I_DA[0]*Demand[t]*P_D[t]+(1-I_DA[0])\
*S_v[t]*P_I[t])-(I_g[t]*P_I[t]*S_g[t]\
+I_d[t]*(FCV[0]*S_d[t]+FCV[1]*Cap)+V_V*S_v[t])
m.Equations(S_v[t]+S_g[t]+S_d[t]-(I_DA[0]*Demand[t]\
+(1-I_DA[0])*S_v[t])==0 for t in range(row))
m.Equations(S_v[t]==(G_w[t]+G_p[t]) for t in range(row))
m.Equations(I_d[t]+I_g[t]-I_DA[0]<=0 for t in range(row))
m.Equations(I_g[t]<(1+S_g[t]-epsilon) for t in range(row))
m.Equations(S_g[t]<(Big_M*I_g[t]) for t in range(row))
#
m.Equations(I_d[t]<(1+S_d[t]-epsilon) for t in range(row))
m.Equations(S_d[t]<(Big_M*I_d[t]) for t in range(row))
m.Maximize(m.sum(c))
m.solve(disp=True)
print('Results')
for t in range(row):
print('I_g{}={}, I_d{}=\
{}'.format(t,str(I_g[t].value),t,str(I_d[t].value)))
print('Objective: ' + str(m.options.objfcnval))
You may also be interested in the m.if3() function such as:
m.Equations([I_g[t]==m.if3(S_g[t],0,1) for t in range(row)])
m.Equations([I_d[t]==m.if3(S_d[t],0,1) for t in range(row)])