teemotheee
teemotheee

Reputation: 53

Dispatch Optimization Problem Using Gekko Cannot Find Solution for Certain Demand Values

I am new to gekko and have been using it to create a plant dispatch scheduler program. The code takes demand values, and schedules the plants while at the same time satisfying the constraints such as dispatch range of each plant, switching states, etc. The intention of the program is to minimize the cost of the objective function which is the total dispatch cost. The program however is not able to solve schedules for all demand values.

The program code is written as follows:

from gekko import GEKKO
import numpy as np
import pandas as pd
# pd.set_option('display.max_columns', None)
# pd.set_option('display.max_rows', None)
m = GEKKO(remote=False) # Initialize gekko
m.options.SOLVER=1  # APOPT is an MINLP solver
m.solver_options = ['minlp_as_nlp 1']
# m.options.IMODE= 3



# demand = [17.23,16.47,15.18,14.42,16.66,19.84,18.24,18.19,26.07,27.34,27.25,26.16,26.01,27.39,28.99,26.83,24.84,27.42,31.23,30.26,29.79,26.73,14.55,20.43]
demand = [17.23]
pmax = [3.75,3.75,3.75,3.25,9.5,7.98,4.6,7.4] #DMCI, POWER ONE, ORMIN, TTR, MHEC #44,480 max
pmin = [2.5,2.5,2.5,2.5,2.48,2.66,0,0]
f_cost = [5000,5000,5000,5000,5000,5000,5000,5000]
meot_cost = [1.8,1.8,1.8,1.8,5,9.8,24.2,3.2]
start_cost = [1000,1000,1000,1000,1000,1000,1000,1000]
v_cost = [2,2,2,2,1,1,1,1]
meots  = [10,10,10,10,10,10,10,10]
n_plants = len(pmax)
hours = len(demand)
# print(len(meots))
print()

# Initialize Variabls
prod = m.Array(m.Var,(hours,n_plants),lb=0, integer=True)
switch_on = m.Array(m.Var,(hours,n_plants),lb=0,ub=1,integer=True)
meot_breach = m.Array(m.Var,(hours,n_plants),lb=0,ub=1,integer=True)   
plant_stat = m.Array(m.Var,(hours,n_plants),lb=0,ub=1,integer=True)   

# Constraints
for i in range(hours):
    #Demand Constraint
    m.Equation(m.sum(prod[i,:]) == demand[i])
for i in range(hours):
    for j in range(n_plants):
        #Plant Dispatch Range Constraint
        m.Equation(prod[i,j] >= pmin[j])
        m.Equation(prod[i,j] <= pmax[j])
        #Plant State constraint
        q = m.if3(prod[i,j],0,1)
        m.Equation(plant_stat[i,j] == q)
        #Switching constraint
        if i == 0:
            m.Equation(switch_on[i,j] == plant_stat[i,j])
        else:
            m.Equation(switch_on[i,j] >= plant_stat[i,j] - plant_stat[i-1,j])
            m.Equation(switch_on[i,j] <= 1 - plant_stat[i-1,j])
        #MEOT Constraint
        running_sum = m.sum(prod[0:i,j])-meots[j]
        x = m.if3(running_sum,0,1)
        m.Equation(meot_breach[i,j] == x)
        
#Objective function
fixed_cost = m.Intermediate(m.sum([plant_stat[t,i]*f_cost[i] for i in range(n_plants) for t in range(hours)]))
start_up_cost = m.Intermediate(m.sum([switch_on[t,i]*start_cost[i] for i in range(n_plants) for t in range(hours)]))
prod_cost = m.Intermediate(m.sum([prod[t,i]*v_cost[i] for i in range(n_plants) for t in range(hours)]))
meot_cost = m.Intermediate(m.sum([prod[t,i]*meot_breach[t,i]*meot_cost[t] for i in range(n_plants) for t in range(hours)]))

total_cost = fixed_cost +\
             start_up_cost +\
             prod_cost +\
             meot_cost
m.Minimize(total_cost)
# m.open_folder()
m.solve(disp=True) # Solve

fixed_cost_value = fixed_cost.value[0]
start_up_cost_value = start_up_cost.value[0]
prod_cost_value = prod_cost.value[0]
meot_cost_value = meot_cost.value[0]

#OUTPUTS
prod = pd.DataFrame(prod)
meot_breach = pd.DataFrame(meot_breach)
plant_stat = pd.DataFrame(plant_stat)
switch_on = pd.DataFrame(switch_on)


print('Production Cost:', prod_cost_value)
print('Startup Cost:', start_up_cost_value)
print('Fixed Cost:', fixed_cost_value)
# print('Meot Cost:', meot_cost_value)
print(prod)
print(plant_stat)
print(meot_breach)
print(switch_on)
print('Objective: ' + str(m.options.objfcnval))

I have tried running the code above and it gave an output dispatch schedules as expected. For a demand value of 17.23 the output is:

Production Cost: 28.57
Startup Cost: 6000.0
Fixed Cost: 30000.0
       0               1       2      3       4               5      6      7
0  [2.5]  [2.5899999996]  [3.75]  [2.5]  [2.48]  [3.4100000004]  [0.0]  [0.0]
       0      1      2      3      4      5      6      7
0  [1.0]  [1.0]  [1.0]  [1.0]  [1.0]  [1.0]  [0.0]  [0.0]
       0      1      2      3      4      5      6      7
0  [0.0]  [0.0]  [0.0]  [0.0]  [0.0]  [0.0]  [0.0]  [0.0]
       0      1      2      3      4      5      6      7
0  [1.0]  [1.0]  [1.0]  [1.0]  [1.0]  [1.0]  [0.0]  [0.0]
Objective: 36028.57

However, when I run certain demand values or the whole 24 hrs, no solution is found. For example if I change the demand value to 14.2, I get the following output:

---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  0.026000000000000023 sec
 Objective      :  36025.14
 Unsuccessful with error code  0
 ---------------------------------------------------
 
 Creating file: infeasibilities.txt
 Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
 @error: Solution Not Found

I would like to understand why it is happening and how to fix it to make it work for all demand values. I am new to gekko and any suggestions would be appreciated.

Upvotes: 3

Views: 92

Answers (1)

John Hedengren
John Hedengren

Reputation: 14376

Set the production >= instead of == to the demand. This allows overproduction when needed. Otherwise, the problem is infeasible.

for i in range(hours):
    #Demand Constraint
    m.Equation(m.sum(prod[i,:]) >= demand[i])

A mixed integer solution is required to use m.if3(). Remove the solver specification to solve as an NLP instead of an MINLP.

#m.solver_options = ['minlp_as_nlp 1']

Here is the full code that solves with variable demand and with a mixed integer solution.

from gekko import GEKKO
import numpy as np
import pandas as pd
# pd.set_option('display.max_columns', None)
# pd.set_option('display.max_rows', None)
m = GEKKO(remote=False) # Initialize gekko
m.options.SOLVER=1  # APOPT is an MINLP solver

# demand = [17.23,16.47,15.18,14.42,16.66,19.84,18.24,18.19,26.07,27.34,27.25,26.16,26.01,27.39,28.99,26.83,24.84,27.42,31.23,30.26,29.79,26.73,14.55,20.43]
demand = [14.2]
pmax = [3.75,3.75,3.75,3.25,9.5,7.98,4.6,7.4] #DMCI, POWER ONE, ORMIN, TTR, MHEC #44,480 max
pmin = [2.5,2.5,2.5,2.5,2.48,2.66,0,0]
f_cost = [5000,5000,5000,5000,5000,5000,5000,5000]
meot_cost = [1.8,1.8,1.8,1.8,5,9.8,24.2,3.2]
start_cost = [1000,1000,1000,1000,1000,1000,1000,1000]
v_cost = [2,2,2,2,1,1,1,1]
meots  = [10,10,10,10,10,10,10,10]
n_plants = len(pmax)
hours = len(demand)
# print(len(meots))

# Initialize Variabls
prod = m.Array(m.Var,(hours,n_plants),lb=0, integer=True)
switch_on = m.Array(m.Var,(hours,n_plants),lb=0,ub=1,integer=True)
meot_breach = m.Array(m.Var,(hours,n_plants),lb=0,ub=1,integer=True)   
plant_stat = m.Array(m.Var,(hours,n_plants),lb=0,ub=1,integer=True)   

# Constraints
for i in range(hours):
    #Demand Constraint
    m.Equation(m.sum(prod[i,:]) >= demand[i])
for i in range(hours):
    for j in range(n_plants):
        #Plant Dispatch Range Constraint
        m.Equation(prod[i,j] >= pmin[j])
        m.Equation(prod[i,j] <= pmax[j])
        #Plant State constraint
        q = m.if3(prod[i,j],0,1)
        m.Equation(plant_stat[i,j] == q)
        #Switching constraint
        if i == 0:
            m.Equation(switch_on[i,j] == plant_stat[i,j])
        else:
            m.Equation(switch_on[i,j] >= plant_stat[i,j] - plant_stat[i-1,j])
            m.Equation(switch_on[i,j] <= 1 - plant_stat[i-1,j])
        #MEOT Constraint
        running_sum = m.sum(prod[0:i,j])-meots[j]
        x = m.if3(running_sum,0,1)
        m.Equation(meot_breach[i,j] == x)
        
#Objective function
fixed_cost = m.Intermediate(m.sum([plant_stat[t,i]*f_cost[i] for i in range(n_plants) for t in range(hours)]))
start_up_cost = m.Intermediate(m.sum([switch_on[t,i]*start_cost[i] for i in range(n_plants) for t in range(hours)]))
prod_cost = m.Intermediate(m.sum([prod[t,i]*v_cost[i] for i in range(n_plants) for t in range(hours)]))
meot_cost = m.Intermediate(m.sum([prod[t,i]*meot_breach[t,i]*meot_cost[t] for i in range(n_plants) for t in range(hours)]))

total_cost = fixed_cost +\
             start_up_cost +\
             prod_cost +\
             meot_cost
m.Minimize(total_cost)
m.solve(disp=True) # Solve

fixed_cost_value = fixed_cost.value[0]
start_up_cost_value = start_up_cost.value[0]
prod_cost_value = prod_cost.value[0]
meot_cost_value = meot_cost.value[0]

#OUTPUTS
prod = pd.DataFrame(prod)
meot_breach = pd.DataFrame(meot_breach)
plant_stat = pd.DataFrame(plant_stat)
switch_on = pd.DataFrame(switch_on)

print('Production Cost:', prod_cost_value)
print('Startup Cost:', start_up_cost_value)
print('Fixed Cost:', fixed_cost_value)
# print('Meot Cost:', meot_cost_value)
print(prod)
print(plant_stat)
print(meot_breach)
print(switch_on)
print('Objective: ' + str(m.options.objfcnval))

Upvotes: 1

Related Questions