Gekko Code for Resource Optimization problem

Question

I am working on resource optimization problem and writing a GEKKO code to solve. Problem statement is as follows: Let assume there are 2 workers and 3 tasks. Each worker gets some rewards for the tasks it picks. Each task can be assigned to only one worker and each worker can only pick one tasks. The objective is to maximize the total rewards obtained by both workers.

m = GEKKO()

# decision variables
x11 = m.Var(value=0, lb=0,ub=1,integer=True) # variable if worker 1 picks task 1
x12 = m.Var(value=0, lb=0,ub=1,integer=True) # variable if worker 1 picks task 2
x13 = m.Var(value=0, lb=0,ub=1,integer=True) # variable if worker 1 picks task 3
x21 = m.Var(value=0, lb=0,ub=1,integer=True) # variable if worker 2 picks task 1
x22 = m.Var(value=0, lb=0,ub=1,integer=True) # variable if worker 2 picks task 2
x23 = m.Var(value=0, lb=0,ub=1,integer=True) # variable if worker 2 picks task 3


# rewards constants
S11 = m.Const(value= 2) # reward if worker 1 picks task 1
S12 = m.Const(value= 5) # reward if worker 1 picks task 2
S13 = m.Const(value= 8) # reward if worker 1 picks task 3
S21 = m.Const(value= 1) # reward if worker 2 picks task 1
S22 = m.Const(value= 5) # reward if worker 2 picks task 2
S23 = m.Const(value= 7) # reward if worker 2 picks task 3

# Equations
m.Equation(x11 + x12 + x13 == 1) # constraints for worker 1
m.Equation(x21 + x22 + x23 == 1) # constraints for worker 2

m.Obj(-(x11*S11 + x12*S12 + x13*S13 + x21*S21 + x22*S22 + x23*S23)) # Objective

m.options.IMODE = 3 # Steady state optimization
m.options.SOLVER = 1        # change solver (1=APOPT,3=IPOPT)
m.solve(disp=True)          # solve locally (remote=False)

print('Results')
print('x11: ' + str(x11.value))
print('x12: ' + str(x12.value))
print('x13: ' + str(x13.value))
print('x21: ' + str(x21.value))
print('x22: ' + str(x22.value))
print('x23: ' + str(x23.value))
print('Objective: ' + str(m.options.objfcnval))

upon running this I am getting output as:

    apm 122.172.80.123_gk_model3 <br><pre> ----------------------------------------------------------------
 APMonitor, Version 1.0.1
 APMonitor Optimization Suite
 ----------------------------------------------------------------




 --------- APM Model Size ------------
 Each time step contains
   Objects      :            0
   Constants    :            6
   Variables    :            6
   Intermediates:            0
   Connections  :            0
   Equations    :            3
   Residuals    :            3
 
 Number of state variables:              6
 Number of total equations: -            2

 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :              4
 
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.00 NLPi:    3 Dpth:    0 Lvs:    0 Obj: -1.50E+01 Gap:  0.00E+00



Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   1.319999998668209E-002 sec
 Objective      :   -15.0000000000000     
 Successful solution

 ---------------------------------------------------
 
Results
x11: [0.0]
x12: [0.0]
x13: [1.0]
x21: [0.0]
x22: [0.0]
x23: [1.0]
Objective: -15.0

The solution looks correct since I have not added the constraint of each tasks to be assigned to be only to 1 worker. But I add following tasks constraints as :

m.Equation(x11 + x21 == 1) # constraint for task1
m.Equation(x12 + x22 == 1) # constraint for task2
m.Equation(x13 + x23 == 1) # constraint for task3

I am getting following error:

    apm 122.172.80.123_gk_model4 <br><pre> ----------------------------------------------------------------
 APMonitor, Version 1.0.1
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :            0
   Constants    :            6

   Variables    :            6

   Intermediates:            0
   Connections  :            0
   Equations    :            6
   Residuals    :            6
 
 Number of state variables:              6
 Number of total equations: -            5

 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :              1
 
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I: -1 Tm:      0.00 NLPi:    2 Dpth:    0 Lvs:    0 Obj:  0.00E+00 Gap:  



    NaN
 Warning: no more possible trial points and no integer solution
 Maximum iterations
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   1.320000001578592E-002 sec
 Objective      :   0.000000000000000E+000
 Unsuccessful with error code            0
 ---------------------------------------------------
 
 Creating file: infeasibilities.txt
 Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
 @error: Solution Not Found
---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
Cell In[40], line 14
     12 m.options.IMODE = 3 # Steady state optimization
     13 m.options.SOLVER = 1        # change solver (1=APOPT,3=IPOPT)
---> 14 m.solve(disp=True)          # solve locally (remote=False)
     16 print('Results')
     17 print('x11: ' + str(x11.value))

File ~/anaconda3/lib/python3.11/site-packages/gekko/gekko.py:2185, in GEKKO.solve(self, disp, debug, GUI, **kwargs)
   2183 #print APM error message and die
   2184 if (debug >= 1) and ('@error' in response):
-> 2185     raise Exception(response)
   2187 #load results
   2188 def byte2str(byte):

Exception:  @error: Solution Not Found

Since this is a very trivial example and the solution does exist. Could you please help in understanding what could be going wrong in the code? and also how do we debug such errors?

Thanks

Answer 1

The second set of constraints should not be equality constraints, but inequality constraints. Reproducing on PuLP,

import pulp

worker_idx = range(1, 3)
task_idx = range(1, 4)

assigns = pulp.LpVariable.matrix(
    name='assign',
    cat=pulp.LpBinary,
    indices=(worker_idx, task_idx),
)

rewards = (
    (2, 5, 8),
    (1, 5, 7),
)

m = pulp.LpProblem(name='resources', sense=pulp.LpMaximize)
m.setObjective(
    pulp.lpSum(
        pulp.lpDot(worker_assigns, worker_rewards)
        for worker_assigns, worker_rewards in zip(assigns, rewards)
    )
)

# Each worker can only pick one task
for worker, worker_assign in zip(worker_idx, assigns):
    m.addConstraint(
        name=f'excl_w{worker}', constraint=pulp.lpSum(worker_assign) == 1,
    )

# Each task can be assigned up to once
for i, task in enumerate(task_idx):
    m.addConstraint(
        name=f'excl_t{task}', constraint=pulp.lpSum(
            worker_assigns[i]
            for worker_assigns in assigns
        ) <= 1,
    )

print(m)
m.solve()
assert m.status == pulp.LpStatusOptimal
for worker, worker_assigns in zip(worker_idx, assigns):
    task, = (
        task
        for task, assign in zip(task_idx, worker_assigns)
        if assign.value() > 0.5
    )
    print(f'Worker {worker}: task {task}')

resources:
MAXIMIZE
2*assign_1_1 + 5*assign_1_2 + 8*assign_1_3 + 1*assign_2_1 + 5*assign_2_2 + 7*assign_2_3 + 0
SUBJECT TO
excl_w1: assign_1_1 + assign_1_2 + assign_1_3 = 1

excl_w2: assign_2_1 + assign_2_2 + assign_2_3 = 1

excl_t1: assign_1_1 + assign_2_1 <= 1

excl_t2: assign_1_2 + assign_2_2 <= 1

excl_t3: assign_1_3 + assign_2_3 <= 1

VARIABLES
0 <= assign_1_1 <= 1 Integer
0 <= assign_1_2 <= 1 Integer
0 <= assign_1_3 <= 1 Integer
0 <= assign_2_1 <= 1 Integer
0 <= assign_2_2 <= 1 Integer
0 <= assign_2_3 <= 1 Integer

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

At line 2 NAME          MODEL
At line 3 ROWS
At line 10 COLUMNS
At line 41 RHS
At line 47 BOUNDS
At line 54 ENDATA
Problem MODEL has 5 rows, 6 columns and 12 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 13 - 0.00 seconds
Cgl0004I processed model has 5 rows, 6 columns (6 integer (6 of which binary)) and 12 elements


Result - Optimal solution found

Objective value:                13.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.00

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.00

Worker 1: task 3
Worker 2: task 2